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Projects  Studied  Under  the  Contract 

During  the  second  half  of  the  contract  year  (Jan.  15,  1972  - 
Jan.  15,  1973),  the  program  continued  the  following  studies: 
speech  analysis  by  linear  prediction,  reconstruction  of  multi¬ 
dimensional  signals  from  projections,  development  of  a  high  speed 
digital  processor  for  speech  synthesis,  and  the  design  of  two- 
dimensional  recursive  digital  filters.  These  projects  are  summar¬ 
ized  in  the  following  rages.  Reprints  of  available  publications 
are  appended . 

The  views  and  conclusions  contained  in  this  document  are  those 
of  the  author,  and  should  not  be  interpreted  as  necessarily  repre¬ 
senting  the  official  policies,  either  expressed  or  implied,  of  the 
Advanced  Research  Projects  Agency  or  the  U.S.  Government. 

1.  Speech  Analysis  by  Linear  Prediction 

Up  to  the  present,  most  of  the  effort  has  been  devoted  to  the 
development  of  an  interactive  speech  analysis  system,  which  is 
implemented  on  the  Fast  Digital  Processor.  The  analysis  procedure 
is  based  on  the  technique  of  linear  predictive  analysis.  In  this 
scheme,  the  combined  effect  of  glottal  source,  vocal  tract,  and 
radiation  losses  is  represented  by  a  single  all-pole  filter.  In 
this  way,  spectral  frames  are  constrained  to  have  a  fixed  number  of 
resonant  peaks,  which  are  located  at  the  positions  of  the  formants 
in  conventional  speech  analysis.  These  spectra  arc  smooth,  yet  they 
retain  the  important  information  about  the  formants.  The  implemented 
system  allows  for  the  rapid  and  interactive  analysis  of  speech 
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samples,  and  the  convenient  storage  of  the  computed  spectral  frames. 

It  is  possible  to  mark  and  edit  the  data  in  flexible  ways,  and  to 
compute  and  attach  auxiliary  calculations  to  the  basic  data.  Any 
of  these  computed  parameters  can  be  easily  displayed  and  examined. 

The  system  will  be  used  to  create  a  very  large  data  base  of 
processed  utterances  for  use  in  speech  recognition  research,  which 
is  the  next  phase  of  this  project.  Disyllabic  utterances  of  the 
form  /h9CiVC2/  will  be  recorded  for  all  possible  C^,  C2  combina¬ 
tions,  and  for  a  wide  range  of  the  vowels  V.  The  goal  will  then 
be  to  describe  the  acoustic  phonetic  parameters  ot  in  these 
utterances.  While  there  are  many  known  contextual  effects  on  the 
phonetic  realization  of  sounds,  it  is  felt  that  these  are  mini¬ 
mized  for  consonants  in  pre-stressed  position.  Furthermore, 
stressed  syllable  nuclei  are  reliable  anchor  points  at  which  to 
initiate  phonemic  recognition.  For  these  reasons,  recognition  of 
pre-stressed  consonants  can  be  expected  to  be  at  least  as  reliable 
as  that  of  consonants  in  other  positions,  and  hence  a  substantial 
effort  toward  their  recognition  is  justified. 

We  expect  that  this  study  will  represent  a  comprehensive 
attack  on  this  problem,  leading  to  the  creation  of  an  exceptionally 
large  data  base,  together  with  a  wide  range  of  techniques  for 
consonant  recognition  and  a  critical  evaluation  of  their  capabilities. 

2.  Reconstruction  of  Multidimensional  Signals  from  Projections 

In  many  applications,  a  set  of  projections  of  an  N-dimensional 
object  onto  (N-l)  dimensions  are  available  from  which  it  is  possible 
to  reconstruct  the  original  object.  X-ray  photographs,  as  obtained 
in  medical  applications  and  the  inspection  of  mechanical  systems, 
represent  two-dimensional  projections  of  the  three-dimensional 
objects  which  have  been  X-rayed.  A  nurnber  of  new  results  for  this 
problem  have  been  obtained  by  formulating  the  reconstruction  problem 
directly  in  digital  signal  processing  terms.  Based  on  thi3  formula* 
tion,  several  algorithms  have  been  developed  which  appear,  based  on 
several  reconstruction  examples,  to  bo  superior  to  previous  algorithms 


in  some  cases.  Most  of  the  reconstructions  performed  have  been 
the  transformation  of  one-dimensional  projections  to  two-dimensional 
photographs.  A  reconstruction  of  a  section  of  leg  bone  from  real 
x-ray  data  has  also  been  performed.  This  work  has  resulted  in  the 
completion  of  an  Sc.D.  thesis  (R.M.  Merserau;  "Digital  Reconstruc¬ 
tion  of  Multidimensional  Signals  from  their  Projections") ,  the 
abstract  of  which  is  attached.  Some  of  the  algorithms  are  also 
summarized  in  the  paper  "The  Digital  Reconstruction  of  Multidimen¬ 
sional  Signals  from  their  Projections"  by  R.M.  Merserau;  Proc.  10th 
Annual  Allerton  Conference  on  Circuit  and  System  Theory,  Oct.  4-5, 
1972,  Monticello,  Ill.,  pp.  326-334. 

Some  preliminary  investigations  into  the  use  of  projections 
for  picture  bandwidth  compression  have  also  been  completed,  which 
have  led  to  promising  results.  If  a  function  can  be  represented  by 
its  projections,  then  perhaps  pictures  can  be  transmitted  or  stored 
by  utilizing  a  set  of  projections  for  these  pictures,  thus  using 
fewer  bits  than  are  required  for  transmitting  the  entire,  picture 
directly.  These  experiments  are  summarized  in  the  doctoral  thesis  by 
Merserau. 

3.  Development  of  a  Digital  Processor  for  Speech  Synthesis 

Detailed  design  of  the  "Black  Box"  processor  has  been  completed. 
There  have  been  two  major  changes  in  design,  and  several  smaller  ones 
as  the  specific  logic  has  been  developed.  First,  a  major  decision 
has  been  made  to  build  the  processor  from  ECL  10k  logic,  rather 
than  74  series  and  Schottky  TTL.  This  change  will  increase  the 
speed  of  the  processor,  but  it  should  also  result  in  greater  system 
reliability.  It  is  our  belief  that  with  proper  design  ; recautions 
as  currently  understood,  ECL  systems  should  be  more  reliable  ar.d 
immune  from  noise  problems  than  TTL  systems.  The  second  major  change 
has  been  an  increase  in  data  word  length  from  18  to  24  bits  to  per¬ 
mit  retention  of  more  significant  figures.  The  new  multiplier  will 
be  16  X  24,  and  should  run  in  less  than  100  nsec.  The  memory  uses 
1024X1  ECL  chips  which  are  just  now  becoming  available,  and  which 
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provide  adequate  storage  at  state-of-the-art  speeds. 

Large  circuit  boards  have  been  ordered,  and  ECL  parts  will  be 
ordered  shortly.  The  only  remaining  design  problems  concern  the 
exact  nature  of  the  PDP-9  interface,  which  must  now  be  modified  to 
allow  for  24-bit  data  words. 

An  internal  document,  "The  Black  Box",  is  attached  which 
describes  the  design  of  the  processor  in  detail. 

4.  Design  of  Two-Dimensional  Recursive  Digitao.  Filters 

Recent  work  has  been  based  on  using  the  one-pro jection  results 
of  Merserau  to  allow  the  inference  of  two-dimensional  structures 
from  their  one-dimensional  projections.  In  this  way,  one-dimensional 
approximation  theorems  can  be  used  to  simplify  the  two-dimensional 
recursive  approximation  problem.  The  theory  for  the  design  of  two- 
dimensional  recursive  filters  whose  magnitude-squared  frequency 
response  approximates  a  desired  function  is  complete,  and  the 
algorithm  has  been  programmed  and  is  currently  being  debugged. 

There  are,  however,  theoretical  problems  concerning  the  realizability 
of  the  designed  filters.  These  problems,  which  must  bo  circumvented 
before  the  design  algorithm  can  be  considered  complete,  are  the  focus 
of  current  effort. 

Additionally,  a  simple  inverse  filtering  program  was  written 
using  the  FDP-Univac  picture  processing  system  developed  last  year. 
Digital  photographs  have  been  blurred  by  a  recursive  lowpass  filter 
(taking  approximately  10  seconds) ,  and  then  reconstructed  exactly 
using  a  non-recursive  high-pass  filter  (taking  approximately  5 
seconds) . 
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ABSTRACT 


Several  algoritlims  for  the  reconstruction  of  multidimensional  signals 
from  their  projections  are  presented.  These  algorithms  can  be  applied 
to  the  problem  of  estimating  the  structure  of  an  unknown  three-dimen¬ 
sional  object  from  its  x-ray  photographs  or  electron  micrographs  taken 
at  different  orientations.  The  reconstruction  problem  is  broken  into 
two  distinct  steps;  first  samples  of  the  Fourier  transform  of  the  un¬ 
known  signal  are  computed  from  a  series  of  digitized  projections,  then 
the  unknown  is  estimated  from  the  samples  of  its  Fourier  transform. 
Reconstructions  are  considered  from  several  sets  of  samples  in  Fourier 
apace.  A  particular  set  of  samples,  the  concentric  squares  raster  is 
developed,  the  reconstructions  from  which  are  superior  to  those  made 
from  the  more  traditional  polar  raster  of  samples  for  bandllmitcd  in¬ 
puts  which  have  a  rectangular  frequency  band.  Furthermore  for  an  im¬ 
portant  class  of  unknowns  exact  reconstructions  can  be  performed  from 
a  concentric  squares  raster  from  a  finite  number  of  projections.  In 
fact  for  this  class  of  unknowns  a  single  projection  is  sufficient.  A 
detailed  treatment  of  the  onc-proj action  reconstruction  problem  is  pre¬ 
sented  and  the  difficulties  associated  with  its  solution  are  explored. 
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THE  DIGITAL  RECONSTRUCTION  OF  MULTIDIMENSIONAL  SIGNALS  FROM  THEIR 
PROJECTIONS. 

RUSSELL  M.  MERSEREAU 
Research  Laboratory  of  Electronics 
Massachusetts  Institute  of  Technology 
Cambridge,,  Massachusetts  02139 

ABSTRACT 

Several  algorithms  for  approximating  a  multidimensional  density 
function  in  terms  of  its  projections  ("x-ray  photographs")  at 
different  orientations  are  presented.  This  is  accomplished  by 
means  of  a  theorem  which  states  that  the  Fourier  transform  of  a 
projection  of  a  function  is  a  slice  (central  section)  of  the 
Fourier  transform  of  that  function.  Special  emphasis  is  given 
to  bandlimited  functions  as  these  are  readily  digitized.  Some 
theorems  are  presented  which  state  that  bandlimited  functions 
of  a  certain  class  can  be  represented  by  a  single  projection. 

INTRODUCTION 

There  are  occasions  when  the  structure  of  a  three-dimensional  object  is 
unknown  and  desired  but  only  two-dimensional  projections  of  that  object  are 
available.  A  transmission  x-ray  photograph  might  represent  such  a  projection. 
A  single  x-ray  photograph,  since  it  is  merely  a  two-dimensional  representa¬ 
tion  of  an  inherently  three-dimensional  structvre,  does  not  contain  all  of 
the  information  that  a  physician  might  want,  since  all  detail  in  a  direction 
normal  to  tha  photographic  plate  has  been  superimposed.  One  solution  to 
this  dilemma  is  to  take  X-ray  photographs  from  different  vantage  points  and 
use  this  set  of  x-rays  to  x '.create  a  three-dimensional  image,  say  in  a 
digital  computer.  In  this  paper  we  will  consider  several  techniques  for  ac¬ 
complishing  this.  In  addition  to  performing  reconstructxon  from  x-rays  [4] 

15] [7] [8),  these  algorithms  are  useful  for  reconstructing  from  electron 
micrographs  [2] [3],  fan-beam  radio  telescope  scans  [1]  and  line  responses 
from  linear  shift- Invariant  optical  systems. 

Because  these  algorithms  will  be  implemented  on  a  digital  comp  iter,  we 
are  constrained  to  work  with  sampled  data.  Our  input  projections  must  be 
sampled  and  our  output  reconstruction  will  consist  of  samples  of  the  unknown 
structure.  As  a  result  this  problem  is  best  considered  as  an  inherently 
digital  one.  This  constrains  us  somewhat  by  limiting  the  class  of  signals 
that  can  be  reconstructed,  but  this  is  not  a  serious  concern  since  most  sig¬ 
nals  that  arise  in  practice  can  be  closely  approximated  by  signals  from  this 
restricted  class.  On  the  other  hand,  by  conntraining  ourselves  to  this  one 
class  of  signals,  we  have  developed  some  extremely  powerful  algorithms  and 
interesting  results. 


THE  GENERAL  RECONSTRUCTION  ALGORITHM 

Up  to  this  point  we  have  assumed  that  projection  functions  ere  simply  the 
mathematical  equivalents  of  x-ray  photographs  or  electron  micrographs.  To 
understand  the  problem  more  fully  we  must  make  that  definition  more  precise. 
Let  us  assume  that  one  unknown  signal  can  be  described  by  an  extinction 
function  f(x,y,a)  and  that  an  x-ray  photograph  is  made  of  the  unknown  by  a 
uniform,  collimated  x-ray  beam  with  intensity  Iq  which  propagates  parallel 
to  the  y-axls.  Then,  ignoring  scattering  effects,  the  observed  intensity 
variation  of  the  x-ray  photograph  can  be  described  by: 

I(x,s)  -  I0  exp{  -/  f(x,y,z)  dy}  (1) 

He  can  then  define  the  projection  function  associated  with  this  orientation 

by 


I 


8 


P0(x,z)  -  dy 


(2) 


The  function  of  f(x,y,z)  might  measure  the  density  of  the  unknown  ee  It  does 
to  some  extent  In  the  case  of  x-rsys  or  It  might  measure  the  extent  of 
staining  of  a  specimen  in  the  case  of  electron  micrographs,  etc. 

To  generalize  eq.  (2),  let  us  assume  that  each  projection  is  taken  normal 
to  the  z-axis  and  that  the  projection  plane,  the  u-z  plane,  meets  the  x-z 
plane  at  an  angle  6.  With  this  notation,  eq.  (2)  corresponds  to  the  projec¬ 
tion  at  angle  3-0°.  We  shall  define  the  projection  at  angle  •  by 


p9(u,z)  ■  I  f (ucosO  -  vsinO,  uain8  ♦  vcos6,  z)dv 


(3) 


If  F(<»x,Uy,Mz)  is  the  Fourier  transform  of  f(x,y,z),  and  if  we  define  a 
plane,  the  w-wz  plane  to  intersect  the  wx-wz  plane  at  an  angle  8,  then  we 
can  further  define  the  slice  of  F(ux,Wy ,wz)  at  angle  8  to  be  F(wcos8,wain8 , 
wz),  i.e.,  a  slice  of  F(ux,Wy,<iiz)  corresponds  to  the  function  evsluated  cm 
a  plane,  which  includes  the  uz  axis,  and  which  is  specified  by  a  single 
angular  parameter.  These  relstlonshlps  are  illustrated  in  Fig.  1. 

Theorem;  (projection/slice  theorem)  The  projection  of  f(x,y,z)  at  angle 
6  to  the  x-z  plane  is  the  Fourier  transform  of  the  slice  of  FOsx.Wy.Ux) 
at  angle  0  to  the  wx-<i>z  plane.  ' 

Proof;  Let  us  define  s  coordinate  system,  the  w,w,wz  coordinate  system, 
which  is  a  rotation  by  0  (radians)  about  the  wz-axis,  of  the  wx,wy,wz 
coordinate  system. 


Then  since 


w  -  w  cosO  +  w  slnO 

«  y 

w  ■  wxsln0  +  WyCOS0 


F(wx,Uy,wz)  -  /  /  /  f (x,y,z)  exp(-j[xwx+yuy+XMzJ)  dxdydz 


(4) 


(5) 


we  can  write 


F(w,w,wz)  -  F(wx,Wy,«x)  -  /  /  /  f(x,y,z)  exp (- j [ xwcos0  ♦  ywslnft 


-  xwsinO  +  ywcosO  +  n*t)]  dxdydz  .  (6) 


Since 


F(wcos0,wsln6,ws)  ■  F(w,0,ws)  -  /  /  /  f(x,y,z)  exp(-j[xwcos0 


by  defining 

u  -  xcosO  +  yslnO 
v  “  -xsine  4-  ycosO 
we  can  observe 
F(wcos0,wsin0,w 


"f*  ywslnO  +  *“ZJ)  dxdydz  (7) 

(8) 


•  •(  * 

z)  ■  F(m,0,wz)  •  /  /  I  /  f (ucos0-vsin0,  u*in0+vcos0 , 

z)dv|exp(-j(u*rfzwz)  )  dudz  (9) 

m  m 

m  /  f  p.(u,z)  *xp(  - j (uw+sw  ) }  duds  (10) 

— «•  —  m 

Q.B.D 
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This  theorem  forms  the  basis  for  our  reconstruction  algorithms.  First 
va  obtain  as  an  input  a  set  of  projections  (say  from  processed  x-rays).  It 
is  most  convenient  if  these  are  all  taken  normal  to  a  single  axis  which  we 
have  defined  as  the  z-axis  in  our  coordinate  system,  although  this  is  not 
necessary  for  the  projection/slice  theorem  can  be  generalized.  Then  we 
transform  these  projections,  by  low  pass  filtering  them,  sampling  them  and 
performing  a  discrete  Fourier  transform  (DFT)  calculation.  These  samples  of 
the  Fourier  transforms  of  the  set  of  projections  provide  us  with  samples  of 
7(wx,<i>y,wz)  by  virtue  of  the  projection/slice  theorem.  From  these  we  can* 
estimate  the  entire  Fourier  space  and  then  by  inverse  Fourier  transforming 
in  three  dimensions,  we  have  an  estimate  of  f(x,y,z),  or,  in  our  terminology, 
a  reconstruction.  How  many  projections  are  needed,  how  they  should  be  sam¬ 
pled  and  how  the  Fourier  space  should  be  estimated  from  samples  of  the 
slices  depend  upon  the  individual  algorithms. 

SPECIFIC  RECONSTRUCTION  ALGORITHMS 

Let  us  assume  that  the  functions  to  be  reconstructed  are  bandlimited. 

.This  assumption  will  guarantee  that  all  of  the  projections  are  bandlimited 
and  it  will  thus  allow  us  o  compute  Fourier  transforms  from  samples  of  the 
projections  with  no  loss  of  information,  a  necessity  for  performing  a  recon¬ 
struction  digitally.  Also  for  our  arguments  let  us  assume  that  we  are 
trying  to  reconstruct  a  two-dimensional  function  (picture)  from  one-dimen¬ 
sional  projections.  This  is  equivalent  to  reconstructing  for  only  a  single 
value  of  z. 

One  approach  to  use  to  obtain  samples  of  the  Fourier  transform  of  an  un¬ 
known  picture  is  to  sample  each  projection  at  the  same  sampling  rate,  that 
rate  being  greater  than  the  maximum  Nyquist  rate  of  each  of  the  projections. 

If  then  each  sequence  of  samples  is  Fourier  transformed  using  a  DFT  algo¬ 
rithm,  the  Fourier  transform  of  the  unknown  picture  f(x,y)  will  be  known  on 
a  polar  raster  of  points.  If  the  projections  were  evenly  spaced  from  0  to 
v,  the  raster  of  points  is  that  shown  in  Fig.  2(a).  From  these  samples  of 
F(mx,w  )  there  are  a  number  of  techniques  that  can  be  used  to  estimate  f(x,y). 
One  of7 the  more  successful  that  we  have  used  is  to  use  linear  Interpolation 
to  estimate  the  values  of  F(ux,Wy)  on  a  Cartesian  (square)  raster  from  the 
polar  samples.  From  the  Cartesian  raster  we  can  perform  an  Inverse  two- 
dimensional  DFT  to  obtain  our  estimate.  Another  approach  that  has  yielded 
good  reconstructions  from  a  polar  raster  is  to  write  the  Inverse  transform 
Integral  in  polar  coordinates, 

1  “  1  - 

F(x,y)  ■  f  I  F(w,9)  exp( J (xwcosfl  +  ywsln9)}w  dwdO  (11) 

0  -* 

and  then  approximate  (11)  by  a  sum,  where  eachr  summand  depends  upon  one  of 
the  polar  samples.  The  latter  technique  performs  well  resolving  detail  but 
does  poorly  on  areas  of  low  Information  content,  such  as  backgrounds  or 
areas  of  nearly  constant  grey  level.  The  interpolation  technique,  on  the 
other  hand,  does  just  the  opposite.  It  resolves  the  flat  areas  but  does 
not  extract  detail  well. 

As  an  alternative  to  sampling  each  projection  at  the  same  sampling  rate 
we  can  vary  the  sampling  rate  with  the  projection  angle.  In  particular, 
asaume  that  f(x,y)  is  bandlimited  within  a  square  in  the  frequency  domain 
such  that  F(wx,wy)  z  0  for  |wx|>W  or  •  Then  the  bandwidth  oi  the 

projection  at  angle  9  is  defined  by 

W  •  - - -  (12) 

max{ |cos9 | ,  |sln6|) 
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How  suppose  we  sssplc  etch  projecdon  At  s  stapling  rate  which  is  propor¬ 
tional  to  its  bandwidth.  If  these  sequences  are  then  DFT'd,  ?(ux,My)  will 
be  known  on  t  concentric  square  raster  such  as  thst  in  Fig.  2(b).  At  first 
glince  such  a  set  of  points  is  not  to  be  favored  over  the  polsr  set,  but 
such  is  not  the  case.  Using  both  the  Interpolation  algorithm  and  the 
integral  approximation  algorithm,  we  get  better  reconstructions  from  the 
concentric  squares  rsstsr  than  from  the  polar  one.  One  resson  for  this  can 
be  seen  from  the  fact  that  the  Cartesian  raster  of  samples  which  we  oust 
have  in  order  to  perform  an  inverse  DPT  lie  along  the  same  horizontal  and 
vertical  lines  as  the  sides  of  the  concentric  squares.  Thus  instesd  of 
the  necessity  of  performing  a  two-dimensional  linear  interpolation,  we  only 
need  to  perform  a  one-dimensional  one.  This  might  be  expected  to  produce 
less  error. 

The  second  advantage  to  a  <  oncentrlc  squares  rsster  becomes  apparent 
when  ve  consider  a  special  cliss  of  signals.  We  have  assumed  thst  f(x,y) 
is  bandllmited.  Let  us  now  assume  thst  it  can  be  represented  in  the  form 


N-l  N-l 


sin  w(x-  — )  alnw(y-  ~) 


*(x.y)  -  l  l  *<y-y“> 

m-0  n-0  w  *w 


W*(x-  =)(y-  =£> 


(13) 


Thus  in  addition  to  requiring  thst  f(x,y)  be  bandllmited,  ve  hsve  required 
that  when  sampled  at  its  Hyqulst  rste  it  have  only  a  finite  number  of  non¬ 
zero  samples.  The  number  N  which  represents  the  width  of  the  square  of 
non- zero  samples  ve  will  refar  to  as  the  order  of  f(x,y).  Bandllmited 
functions  of  finite  order  are  completely  specified  by  their  DFT's  and  their 
Foutlcr  transforms  arc  two-dimensional  polynomials  of  degree  N-l  in  each 
variable.  Because  of  the  fact  that  a  one-dimensional  polynomial  of  degree 
d  is  completely  specified  by  d+1  independent  samples,  it  can  be  proven  that 
a  bandllmited  function  of  order  N  can  be  reconstructed  exactly  from  N+l 
projections  (although  not  necesssrlly  for  all  sets  of  N+l  projections). 

Thus  for  this  class  of  functions,  concentric  square  projections  as  these 
projections  shall  be  called,  assume  theoretical  importance. 

In  Fig.  3  is  shown  a  reconstruction  of  a  cross  section  of  a  leg  bone 
from  36  concentric  squares  projections.  The  projections  were  secdons  of 
x-rsys  which  were  taken  at  5°  intervale  which  are  scanned  logarithmically. 
Each  section  of  esch  x-rsy  consisted  of  236  samples  and  the  reconstruction 
is  plotted  on  a  236x236  Cartesian  raster. 


RECONSTRUCTING  FROM  ONE  PROJECTION 

In  the  previous  section  we  noticed  thst  bandllmited  functions  of  finite 
order  could  be  reconstructed  exactly  from  N+l  (concentric  squares)  projec¬ 
tions,  where  N  was  the  order  of  the  function.  *  Is  this  the  minimum  number 
of  projections  needed  to  reconstruct  thess  functions?  The  answer  to  this 
question  is  no.  In  fact,  functions  of  this  form  can  be  reconstructed  from 
a  single  projection. 

Theorem  (one  projection  theorem)  A  bandllmited  function  f(x,y)  of 
order  N  in  each  variable  can  bs  reconstructed  from  the  concentric 
squares  projection  at  8  •  tan”l(l/N). 

Proof:  Consider  the  slice  at  8  •  tan"^(l/H). 

F(«cos8,wsln8)  -  F  (  - ■  ,  - — - )  (16) 

All  +  1  /N2  +  1 

where 
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I 


N-l  N-l 


F(V",>  -  Jo  Jo  <<f-  ’  f)  i(cV™y)])  bww<“x’“y) 


(15) 


where 


b..,(iu  ,u  )  “ 

WW  x*  y 


fl  .  I“xl  IW,  |u»y|  <  W 
0  ,  otherwise 


Evaluating  (14), 


nu 


F(  ,  -2-  )  -  Y  Y  £(—  *  22-)  exp{-j  (- 

'  4P+I  m-0  n-0  W  W  Wr^ 


[N»+n])>  (16) 


|b)|  <.{} 


If  we  define  g(Nm+n)  ■  f(— ,  21)  then  (16)  becomes 

W  U 

F(  ■■  ■“—  ,  ~ —  )  -  l  g(p)  exp(-J  •  *“P~ 

M  /nJ+T  p-o  w/nmT 


),  |M|  /H2TI 


otherwise 


(17) 


Thus  over  the  region  of  interest,  the  slice  at  6  ■  tan"^(l/H)  (which  can 
be  obtained  from  one  projection)  is  a  one-dimensional  polynomial  of  degree 
N^-l  in  the  variable  exp{-j  (xui/w/n^+1)  ) ,  and  the  coefficients  of  that 
polynomial  are  simply  the  picture  samples  arranged  column  by  column.  Thus 
knowledge  of  sample  values  specifies  the  picture  samples  and  by  (13) 
this  specifies  the  unknown  picture  function. 

Q .E. D » 


This  theorem  also  has  implications  in  two-dimensional  filter  design  for 
it  implies  that  the  frequency  response  of  a  two-dimensional  non-recursive 
digital  filter  is  specified  by  its  values  along  a  single  radial  line, 
although  how  this  fact  might  be  utilized  in  filter  design  is  still  not 
clearly  understood. 

Despite  the  beauty  of  the  one-projection  theorem,  it  is  not  particularly 
useful  as  a  reconstruction  technique  on  a  finite  precision  machine,  because 
for  values  of  N  larger  than  4  or  5,  solving  eq.  (17)  for  {f(mir/W,  nx/VO) 
requires  the  inversion  of  a  large  nearly  singular  matrix  and  any  errors 
made  in  obtaining  the  slice  samples,  which  are  bound  to  occur,  are  amplified 
enormously.  In  this  respect,  applying  the  one  projection  theorem  is  much 
like  trying  to  apply  analytic  continuation  to  an  unknown  function  all  of 
whose  derivations  are  known  at  a  single  point. 

From  a  theoretical  viewpoint,  it  is  interesting  to  ask  if  there  are  any 
other  of  theae  critical  slices.  Clearly  not  all  slices  will  work;  for  ex¬ 
ample,  knowledge  of  sample  values  along  Just  the  ux  or  uiy  axis  is  generally 
not  sufficient.  In  fact,  there  are  an  infinite  number  or  these  critical 
slices.  A  necessary  and  sufficient  condition  for  a  slice  with  a  rational 
slope  to  be  a  critical  slice  is  given  in  the  following  theorem  which  will 
be  stated  without  proof. 


Theorem:  If  a  slice  has  an  angle  tan~^(A/B)  where  A  and  B  arc 
integers  with  no  common  factor,  then  a  necessary  and  sufficient 
condition  for  this  slice  to  be  sufficient  for  reconstruction  of 
a  bandlimited  function  of  order  H  is  that  the  equation 

Bm  +  An  -  Bm*  +  An*  ,  0  <_  m,n,m'  ,n'  <_  N-l 
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possess  only  the  trivial  solution  m  -  m* 

n  ■  n’ . 

In  particular,  if  N  is  a  power  of  two  and  A  is  an  odd  integer,  the 
slice  with  slope  tan"l(A/N)  is  critical. 

Using  techniques  completely  analogous  to  that  by  which  we  derived  the 
one  projection  theorem,  we  can  derive  two  projection  theorems,  four  pro¬ 
jection  theorems,  etc.  As  a  rule  of  thumb,  the  greater  the  number  of  pro¬ 
jections  which  we  care  to  use,  the  easier  it  is  to  obtain  a  solution  with 
real  data.  It  generally  takes  about  N/2  projections  for  the  sensitivity  to 
be  reduced  enough  to  use  currently  available  machines. 

SUMMARY 

A  number  of  algorithms  have  been  presented  for  reconstructing  multi¬ 
dimensional  signals  from  their  projections,  such  as  x-rays.  In  addition  to 
their  uses  in  performing  reconstructions,  projection  functions  are  poten¬ 
tially  useful  for  characterizing  multidimensional  signals  for  purposes  of 
pattern  recognition  or  bandwidth  compression  for  signal  transmission.  Uses 
"for  projections  in  these  areas  has  not  been  explored. 
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THE  BLACK  BOX 
being  a 

small,  fast,  inexpensive 
digital  processor 
desigried  mainly  for 
speech  synthesis 
but  nicely  suited 

t 

for 

myriad  other  tasks 

J.  Allen 
F.X.  Carroll 
E.  Jensen 


This  document  contains  a  detailed  description 
of  the  proposed  Black  Box.  Readers  are  asked 
to  comment  freely  and  quickly  so  we  can  proceed 
with  construction. 
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The  Black  Box 


1 •  Introduction.  In  this  paper,  we  will  describe  a  new  computer 
having  several  unusual  design  features.  The  original  motivation 
for  this  design  was  the  need  for  a  real-time  all-digital  speech 
synthesizer.  Since  the  vocal  tract  model  to  be  simulated  is 
complicated  (see  Fig.  1,  due  to  D.H.  Klatt) ,  it  was  necessary 
to  adopt  several  features  which  optimize  these  calculations  in 
time.  The  heart  of  the  computer  is  a  very  fast  multiplier 
(18  X  18  in  about  100  nsec) .  No  instruction  overlap  is  used, 
but  instructions  have  a  three-address  format,  so  that  (for  example) 

A-B  -*•  C  is  done  in  one  instruction,  including  the  two  loads  and 
one  deposit.  Each  instruction  contains  an  op-code  plus  these 
three  addresses,  as  well  as  a  jump  address  for  the  next  instruc¬ 
tion,  so  that  no  program  counter  is  needed.  There  is  a  separate  instruction 
memory  of  44-bit  width,  and  two  data  memories,  each  of  18-bit 
width.  It  should  be  noted  that  no  special  registers  are  provided. 

There  is  no  accumulator,  no  program  counter,  and  the  machine 
status  word,  as  well  as  the  direct  memory  access  word  count  and 
addresses  are  kept  in  memory,  so  that  they  can  always  be  inspected 
by  the  program.  Very  little  timing  is  needed  internal  to  an 
instruction,  since  the  multiplier  is  combinatorial,  and  shifting 
is  accomplished  via  multiplexing.  This  results  in  a  fairly  simple 
control  structure,  most  of  the  complication  arising  from  I/O 
considerations.  The  computer  is  designed  to  be  interfaced  to  a 
host  PDP-9  rrachine,  and  will  probably  not  be  used  in  a  stand-alone 
mode,  although  this  is  possible.  Direct  memory  access  transfers 
are  possible  in  either  direction  between  the  PDP-9  and  any  memory 
location  in  the  black  box  at  a  1  megacycle  rate.  In  addition, 
transfers  to  and  from  the  PDP-9  accumulator  and  any  black  box 
location  are  possible.  PDP-9  IOT  instructions  can  be  used  to 
control  various  control  bits  in  the  black  box.  Finally,  part  of 
the  data  memory  is  paralleled,  so  that  it  is  possible  to  compute 
using  one  parameter  set,  while  a  new  set  is  being  loaded  (trans¬ 
parently  to  the  black  box  program)  from  the  PDP-9.  In  this  paper, 
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we  present  several  examples  to  show  the  utility  of  this  device 
for  other  specialized  tasks,  including  display  processing, 
floating  point  calculations,  and  FFT's.  The  basic  design  goal 
has  been  to  achieve  a  powerful,  fast,  yet  inexpensive  processor, 
but  with  little  consideration  given  to  ease  of  programming. 

2.  Architecture.  Figure  2  shows  a  rough  block  diagram  of  the 
computer.  Three  components  can  be  recognized. 

^  Instruction  processor?  Includes  a  256  X  44-bit 
instruction  memoty,  loadable  directly  from  the  host  machine; 
an  instruction  register;  instruction  addressing;  skip  logic; 
and  instruction  decoding. 


Data 

a) 


processor:  Includes 

Two  data  memories,  X  +  Y,  each  256  X  18,  and 
each  having  some  of  its  locations  paralled  by 
additional  memory  locations. 

Three  processing  units: 

1.  Multiply  unit,  performs  X.Y+z 

2.  Arithmetic-Logical  Unit  (ALU),  performs 
adds,  subtracts,  and  logical  operations. 

3.  BIT  test,  performs  skip  on  a  selected 
bit  of  a  given  word  and  provides  for 
modification  of  that  bit. 

a  Select:  Selects  desired  output  from  pro¬ 
cessing  unit,  with  capability  to  incorporate  shifts. 


3)  Input  Output  Processor:  provides  programmed  and  DMA 
transfers  with  the  host  computer. 

^-Ure  1  shows  a  more  detailed  block  diagram,  which  is 
intended  to  relate  to  the  further  detailed  discussion. 


3*  Il>-g.truction  processor:  Instructions  are  44  bits  long,  and 
cannot  be  modified  within  the  black  box  except  that  the  effective 
address  of  all  address  fields  can  be  modified  by  index  registers, 
or  an  instruction  may  be  skipped.  The  host  computer,  however, 
has  access  to  the  256-location  instruction  memory  so  that 
instructions  may  be  updated  at  any  time  19 
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by  input  transfers  into  the  black  box.  There  is  no  program  counter# 
so  that  each  instruction  contains  a  jump  address  to  the  next  instruc¬ 
tion.  This  jump  address  can  be  modified  by  indexing  or  by  ORinq  into 
its  LSB  the  OR  of  the  following  3  conditions  t 

1)  Data  Select  output  *  zero 

2)  "  M  "  negative 

3)  one  bit  selected  by  BIT  instruction. 


There  are  five  classes  of  instructions, 
block  of  features: 

4-  .  Z  .  5  7-8  8 


T~ 


All  five  have  a  common 

Q,  vT  0U.A  S 'V ^ 

*  8 


Y\t 


9a 


2L 


Y 

at 


The  only  variations  in  this  format  over  the  entire  group  of 
instructions  are  the  length  of  the  X-address  field  (8  or  9  bits) 
and  the  specifics  of  the  "Y  or  special'*  field#  which  is  either  the 
Y-address  (in  double-source  instructions)  or  the  specifier  of 
actions  peculiar  to  certain  single-source  instructions.  We  first 
describe  the  common  features  (i.e.  everything  but  the  "Y  or  special" 
field)  : 


1)  Index  control: 
the  X  memory 

X-Address 

U8 

128 

138 

148 


There  are  4 
Register 


index  registers  located  in 
Action 

modifies  Y-address 
"  X-  " 

M  D-  " 

"  Jump  " 


The  4  index  bits  control  the  ORing  of  these  registers  into 
their  respective  address  registers.  The  use  of  ORing#  as 
opposed  to  adding,  simplifies  the  logic  and  increases  the 
speed  of  indexing  at  very  little  cost  in  programming  conveni¬ 
ence.  In  single-argument  and  BIT  instructions,  the  Y  index 
bit  in  ignored. 

20 
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2)  Skip  bits:  There  are  always  two  possible  skips  associated 
with  each  instruction. 

N  skip  on  data  select  output  negative 
Z  skip  "  "  "  "  ■  zero 

Placing  l's  in  these  fields  enables  the  skips,  which  cause 
a  1  to  be  ORed  into  the  jump  address  LSB  if  the  designated 
skip  condition  is  met. 

3)  Op  Code:  The  op  code  is  5  bits  long.  Detailed  explan¬ 
ations  are  given  below. 

4)  X-address .  In  double-source  instructions,  this  is  an 
8-bit  address  in  the  X-memory.  In  single-source  instruc¬ 
tions,  the  9  bits  designate  a  location  in  X  or  Y  memory. 

If  the  MSB  is  0,  the  X-memory  is  addressed,  and  if  it  is  a 
1,  the  Y-memory  is  addressed. 

5)  D-address ;  This  is  the  9  bit  destination  address, 
i.e.  where  the  result  of  an  instruction  is  stored.  When 
the  memory  is  in  the  serial  mode,  the  result  can  be  stored 
in  any  X-  or  Y-location,  but  when  the  memory  is  in  parallel 
mode,  the  result  is  stored  in  the  analagous  locations 5 
specified  by  the  8  least  significant  bits,  in  both  memories. 

t 

This  last  situation  is  violated  in  input  transfers,  discussed 
below. 

6)  Jump  Address:  8-bit  address  of  next  instruction. 

This  block  of  features  is  augmented  in  different  ways  to  form 
the  five  instruction  classes: 

I)  Double-Source  Instructions:  The  X  and  Y  fields  are  both 
8  bits  long,  and  each  specify  a  source  location  in  their 
respective  memories.  See  Figure  4.  The  result  is  stored 

in  D,  and  the  next  instruction  taken  from  J. 
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II)  Mi^ltijuly :  X  and  Y  are  8-bit  source  addresses,  and  a 
third  source  is  always  implied,  namely  the  Z-rogister,  which 
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is  location  lOg  in  the  X-memory.  See  Figure  5.  Note  that 
Z  is  always  cleared  before  a  deposit  is  made,  so  that  it 
will  either  be  properly  updated  when  D  -  Z,  or  sec  to  zero 
to  proviue  only  the  X«Y  function.  The  op-code  contains 
*  3-bit  scaling  field,  and  the  8  possible  options  are  shown 
in  Figure  5. 

III)  Index  Register  Modification  Instructions;  These 
look  like  an  ADD  instruction  in  that  the  contents  of  the 
two  source  addresses  are  added  and  stored  in  the  destina¬ 
tion  address,  except  that  the  result  is  also  stored  in  the 
index  register  referred  to  by  the  instruction.  See  Figure  6. 
There  are  4  index  instructions,  one  for  each  index  register. 

IV)  BIT  instruction:  The  X-field  is  9  bits  long  so  that 
any  X  or  Y  location  can  be  used  as  source.  The  remaining 
(special)  field  of  7  bits  is  divided  into  a  5-bit  select 
field  and  a  2-bit  modification  field.  See  Figure  7.  The 
select  field  specifies  which  of  the  18  bits  of  a  data  word 
is  to  be  examined.  If  this  bit  is  a  1,  it  is  ORed  into  the 
LSB  of  the  jump  address.  The  modification  field  can  then 
affect  this  bit  in  any  of  4  ways,  shown  in  the  figure. 
Additionally,  after  this  modification  is  made,  the  usual 
skip  tests  can  be  performed  on  the  resultant  data  word. 

Single-Source  Instructions;  Here,  the  X-address  field 
is  9  bits  long,  and  the  remaining  7  bits  is  used  to  specify 
shifts  and  rotates  as  shown  in  Fic ire  8.  Thus  both  source 
and  destination  can  be  any  address  in  X  or  Y. 

The  instruction  processor  has  the  following  paths  to  the 
rest  of  the  machine: 

a)  Index  bits:  cause  index  registers  to  be  ORed 

into  the  corresponding  address  registers 
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b)  Skip  bits:  cause  outputs  from  skip  tests  to  be 

ORed  into  the  jump  adr.ress 

c)  Op  code  including  "special"  field  when  present: 

goes  to  instruction  decoder,  then  to 
instruction  execution  control  in  the 
data  processor. 

d)  Source  and  destination  addresses:  routed  to 

respective  data  memory  address  registers. 

e)  Jump  address:  routed  to  instruction  address 

register. 


4.  Addressing  Figure  9  shows  the  complete  address  space,  both  as 
seen  internally  and  by  the  host  computer.  Note  that 

1)  To  the  host  machine 

a)  the  X  t  Y  memories  are  one  512-register  memory, 
locations  0-777 

b)  each  instruction  memory  location  is  mapped  onto 
4  host  computer  18-bit  words,  so  that  4X256  -  1024 
host  memory  words  are  needed  to  represent  tho  instruc¬ 
tion  memory.  These  are  locations  2000-3777  of  the 
overall  address  space. 

2)  The  top  408  locations  of  the  Y  memory  can  be  switched 
between  two  separate  physical  memories.  More  on  this  below. 

3)  The  bottom  15g  locations  of  the  X  memory  aro  special 
registers,  which  can  be  accessed  directly  in  addition  to 
the  normal  X-addrcssing. 
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5.  Data  Processor:  As  mentioned  above,  there  are  3  units  in  the 
data  processor: 

a)  X  ♦  Y  data  memories:  These  are  256  X  18-bit  memories, 
which  can  run  in  one  of  2  inodes  selected  by  a  control  bit 
in  Word  0  of  the  X-memory. 

1)  Serial :  This  mode  logically  places  the  X  and  Y 
memories  end-to-end.  X  is  the  bottom  256  locat.'.ons, 
and  Y  the  top  256  locations.  Thus,  in  this  mode,  all 
deposits  and  single-source-instrv.tion  loads  have  access 
to  all  512  locations.  Additionally,  the  I/O  instruc¬ 
tion  (which  is  a  single-source  instruction)  can  access 
all  512  locations  in  this  mode.  Of  course,  double¬ 
source  instructions  access  both  X  ♦  Y  separately  and 
simultaneously. 


2)  Parallel :  Source  loads  are  unaffected  by  this 
mode,  but  deposits  qo  to  analogous  locations  in  both 
X  ♦  Y  memories.  In  the  following  diagram  the  case 
of  the  shadow  memory  in  parallel  mode  is  shown.  When 
ti.e  shadow  mode  and  parallel  mode  are  both  •"Elected, 
input  and  output  transfers  within  the  shadow  memory 
range  go  to/froro  that  memory  which  is  selected 


fVfe&OVl) 


lw«*utvU»^  \ 


as  the  I/O  shadow  memory  at  tne  time  of  the  transfer, 
and  no_- whore  clno.  It  in  as  thouqh  the  other  j>art  of 
the  parallel  transfer  wont  to  a  (non-existent)  X 
shadow.  Thus,  in  qencral,  oven  in  the  parallel  mode, 
the  contents  of  X  i  f  within  the  shadow  ranee  will  not 
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agree,  and  this  is  desirable.  Switching  of  the  shadow 
memories  only  effects  which  of  the  two  is  currently 
local  to  the  black  box  Y  memory,  and  which  is  used 
for  I/O.  (See  below) 

The  X  memory  has  its  bottom  15g  registers  used  for  special 
purposes,  as  shown  in  Figure  10.  These  registers  can  be  used  just 
like  any  other  memory  register,  but  in  addition,  registers  4-14 
have  special  direct  output  lines,  and  registers  0-3  have  direct 
input  and  output  lines  to  the  host  computer.  Additionally,  loca¬ 
tions  1,  2,  3,  and  7  are  counters.  All  of  these  special  registers 
*re  treated  in  detail  below.  Word  O  is  the  control  register,  and 
is  also  given  special  treatment. 

The  Y  memory  has  its  top  321Q  registers  duplicated  by  a 
"shadow"  memory.  (See  Figure  3)  When  the  shadow  mode  is  off,  the 
Y  memory  is  a  straightforward  256  X  18  memory.  But  when  the  shadow 
mode  is  on,  a  bit  in  the  control  register  selects  which  of  the  two 
memories  will  accommodate  local  loads  and  deposits  within  this 
address  range.  I/O  with  the  host  computer,  however,  will  utilize 
that  memory  not  so  selected,  and  this  is  done  automatically.  This 
feature  was  provided  to  allow  for  a  parameter  memory,  one  part  of 
which  could  be  used  for  local  computing,  while  the  other  is  being 
updated  by  the  host  computer. 

b)  Processing  Units:  Within  the  data  processor,  there 

are  three  processing  units,  specialized  by  function. 

1)  Multiply  As  described  in  Figure  5,  this  unit 
performs  X*Y+Z.  Z  is  supplied  by  register  10  of  the 
X  memory,  and  is  always  used  in  the  calculation.  An 
18X18  array  multiplier  perform  full  2's  complement 
multiplication,  together  with  the  Z  add,  in  about 
100  nsec.  All  scalings  are  done  in  the  data  select 
unit. 
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2>  Arithmetic-Logical  Unit  (ALU)  This  is  a  processor 
accepting  X  t  Y  inputs  and  performing  adds,  subtracts, 
and  logical  operations,  it  is  realized  in  74S181's 
and  is  used  for  all  instructions  except  multiply. 

— *t..tC3t  (BIT)  This  unit  selects  one  among  the 
18  bits  of  a  data  word,  and  connects  it  to  the  jump 
address  skip  logic.  Thus,  the  selected  bit  is  Ored 
into  the  LSB  of  the  jump  address.  The  selected  bit 
can  also  be  modified,  as  shown  in  Figure  7,  and  this 
is  accomplished  in  the  ALU. 

c)  Data  Select;  This  unit  is  a  largo  multiplexor,  capable 
of  gating  any  one  of  18-bit  words  through  the  unit. 

These  words  include  input  transfers,  multiplier  scalings, 
and  ALU  shifts  and  rotates.  Note  that  all  shifting  is  done 

by  multiplexing,  so  that  no  shift  register  or  counter  is 
required. 


Skip  tests  for  negative  and  zero  arc  made  at  the  output  of 
the  data  selector,  and  all  result.-,  are  held  in  a  latch, 

while  the  destination  address  is  set  up,  before  they  ore 
stored. 


6-  ineut/output  Processor,  This  processor  handles  all  I/O  trans¬ 
fer.  with  the  host  computer.  As  shown  in  rigurc  9,  the  host 
machine  has  access  to  all  memory  locations  in  the  black  box.  The 
/  processor  accomplishes  each  transfer  by  executing  a  one-instru 
on  nterrupt  (see  bolow) ,  which  is  transparent  to  the  currently 
running  program.  This  is  the  only  interrupt  facility  in  the  black 

box.  Transfers  may  take  place  in  either  direction,  and  are  of  two 
types. 

a)  AC-trnn«forjs.  When  the  PDP-9  is  the  host,  transfers  to 
and  fro.n  its  AC  c.in  be  mado  under  host  I/O  control.  Thus 
these  transfers  can  only  bo  originated  from  the  host  machine 
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b)  DMA  transfers.  DMA  transfers  (up  to  a  5  megacycle  rate) 
can  be  Initiated  by  either  machine.  The  required  special 
registers  are  in  the  X-memory. 

X-Location  Use 

1  Host  DMA  address 

2  Local  DMA  address 

3  Word  count 

Note  that  all  three  of  these  registers  are  counters. 

The  details  of  these  transfers  are  as  follows:  (IOT  refers  to  host 
I/O  control  instructions) . 

a)  AC  transfers  (LOC  refers  to  an  arbitrary  black  box  data 
location) 

1.  PDP-9  AC  — *  BB  LOC 

a)  IOT  puts  BB  address  into  BB  X 2,  via  input 
buffer  register  (IBR) 

b)  IOT  puts  data  in  IBR,  and  causes  interrupt  to 
transfer  data  to  LOC 

2.  BB  LOC  -4  PDP-9  AC 

a)  IOT  puts  BB  address  into  BB  X2,  and  then 
causes  C(X2)  to  be  loaded  into  Output  Buffer 
Register  (OBR) 

b)  IOT  calls  for  data  from  OBR  —4  AC 

b)  DMA  transfers 

1.  PDP-9  DMA  to  BB 

a)  IOT  host  starting  address  to  X^ 

b)  IOT  black  box  starting  address  to  X2 

c)  IOT  word  count  to  X^ 

d)  IOT  to  initiate  transfer:  direction  of  transfer 
and  shadow  mode  control  arc  also  specified. 

2.  BB  to  PDP-9  DMA 

This  is  similar  to  the  above,  except  for  direction 
of  transfer. 

3.  DMA  transfers  can  also  be  initiated  directly  within 
the  black  box,  since  all  those  special  registers, 
and  the  DMA  control  flags  are  within  the  X-memory. 
Obviously,  these  registers  can  be  inspected  at  any 
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The  I/O  instruction,  executed  during  the  I/O  interrupt,  is  a 
somewhat  specialized  MOVE  (See  Figure  8).  Thus  it  has  9-bit 
source  and  destination  addresses. 

a)  Input  LOC  is  in  the  DMA  BB  Address  register  (X  memory 
word  2);  data  select  from  Input  Buffer  register  is  enabled, 
and  a  MOVE  takes  place  from  LOC  to  LOC.  Thus  the  initial 
contents  of  LOC  are  placed  on  the  X-bus,  but  the  outputs 
of  the  data  processors  are  not  gated  through  data  select. 

In  this  way,  the  initial  contents  of  LOC  are  lost,  and  the 
Input  transfer  only  uses  the  "store"  part  of  the  MOVE 
instruction. 

k)  output  LOC  is  in  the  DMA  BB  Address  register.  A  MOVE 
from  LOC  to  LOC  is  executed,  and  the  output  of  the  data  selector 
(the  contents  of  LOC)  is  also  placed  in  the  Output  Buffer 
Register  (OBR) . 

7*  Special  Registers.  As  previously  mentioned,  the  bottom  15g 
locations  of  the  X  data  memory  are  special,  in  that  they  are  more 
chan  plain  memory  locations.  Some  have  special  input  and/or  output 
lines,  and  may  even  be  counters.  The  several  subsets  of  these 
registers  require  special  discussion. 

a)  iilMg*  These  are  the  4  index  registers.  The  index 
instructions  (Figure  6)  store  into  them  as  well  as  the 
designated  destination  address.  The  outputs  of  these 
registers  are  ORed  into  their  respective  address  fields  if 
gated  by  l's  in  the  respective  index  control  bits  of  the 
instruction  being  executed.  These  registers  are  not  counters. 
All  index  incrementing  must  be  performed  explicitly. 

b)  10g .  This  is  the  Z-register .  Its  output  lines  are 
permanently  connected  to  the  multiplier,  so  that  a  multiply 
instruction  always  adds  the  contents  of  Z  to  X*Y  to  form 
the  output  of  the  multiplier. 

c)  ^~7  These  are  the  Clock  and  Count  registers, 
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respectively.  Normally  the  clock  interval  is  in  the  clock 
register.  This  count  is  automatically  jammed  into  the  count 
register  when  the  latter  goes  to  0  from  -1.  This  condition 
is  equivalent  to  the  transition  to  0  being  coupled  with  a 
carry  out  of  the  MSB.  Thus  no  jamming  takes  place  when  0  is 
deposited  into  the  count  register.  Each  time  the  desired  inter¬ 
val  is  jammed  into  the  count  register,  a  pulse  is  generated 
for  external  use.  The  count  register  is  also  constructed 
so  that  it  cannot  count  beyond  0.  If  a  periodic  output 
pulse  is  desired,  merely  set  the  desired  period  in  the  clock 
register,  and  the  count  register  will  generate  output  pulses 
at  the  period  interval  repeatedly.  An  external  clock  supplies 
the  pulse  train  for  the  count  register.  Any  existing  .count 
can  be  overridden  at  any  time  by  simple  depositing  into  the 
count  register.  This  scheme,  of  course,  also  provides  for 
changes  from  one  period  to  another,  by  merely  updating  the 
clock  register  at  the  appropriate  time,  i.e.  before  the 
current  period  expires.  Also,  it  is  clear  that  no  special 

instructions  are  needed  for  these  operations. 

•  *. 

d)  4-5  These  are  the  Output  Registers.  Their  outputs 
are  brought  directly  (buffered,  of  course)  to  rear-panel 
connectors.  The  intent  is  to  connect  D/A  converts  to  these 
outputs . 

e)  1-3  These  are  the  DMA  address  and  word  registers. 

Location  1  contains  the  host  computer  current  address, 
location  2  the  local  (black  box)  current  address,  and  loca¬ 
tion  3  the  present  word  count.  This  last  location  generates 
a  pulse  on  the  transition  from  -1  to  0  for  use  in  supplying 
an  interrupt  to  the  host  computer.  All  three  registers 

are  counters. 
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f)  0  This  is  the  Control  Register.  Its  individual  bits 
are  employed  as  follows:  (no  significance  is  attached  to 
the  order.) 

1.  (1  bit)  Run/Halt.  DMA  will  still  transfer  during 
Halt. 

2.  (1)  Link.  This  flag  is  set  on  carry  out  of  the 
MSB  during  ADD  or  SUB.  There  are  corresponding 
instructions  which  provide  ADD  with  Link  and  SUB 
with  Link. 

3.  (1)  Overflow.  This  flag  is  set  on  overflow  during 
ADD  and  SUB  (and  possibly  on  some  multiplier  outputs) . 
It  can  be  reset  only  by  program  control  (using  BIT). 

4*  (2)  Shadow.  One  bit  turns  the  shadow  mode  on,  so 

that  I/O  transfers  use  the  memory  not  selecte'd,  and 
the  other  bit  selects  which  physical  shadow  memory 
is  used  by  the  BB  for  local  computation. 

5*  (1)  Series/Parallel .  In  the  serial  mode,  deposits 

are  made  to  one  of  512  locations  in  X  or  Y,  whereas 
in  parallel  mode  deposits  go  to  two  locations,  one 
in  each  of  X  and  Y.  (Note  comments-  above  concern¬ 
ing  I/O  in  parallel  mode  within  the  shadow  address 
range) . 

6.  (3)  DMA  There  are  3  control  bits 

a)  Start.  Resets  when  transfer  is  done,  and 

initiates  interrupt  (when  enabled) 
to  the  host  computer  when  done. 

b)  Direction. 

Into  or  out  of  the  Black  Box. 

c)  Interrupt  Enable. 

Enable  flag  allows  Black  Box  to 
interrupt  host  computer  when  DMA 
transfer  is  done. 

7.  (1)  Clock:  This  is  an  enable  flag,  which  allows  the 
counter-register  — »  0  condition  to  interrupt  the 
host  computer. 
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8.  (1)  Program  Interrupt.  This  is  a  general  flag 

which  can  be  used  to  interrupt  the  host  computer 
for  any  reason. 

(1)  Program  Interrupt  Enable.  Provides  overall 
control  of  all  interrupts  from  the  Black  Box  to 
the  host  computer. 

Items  1-9  above  account  for  12  bits,  so  there  are  still  6  bits  left 
which  can  be  used  for  flags  or  other  purposes. 

8.  Front  panel:  The  front  panel  will  have  two  sets  of  switches 
and  corresponding  lights.  One  set  is  18  bits  long  f?r  data  or 
instructions.  The  other  set  is  12  bits  long  and  specifies  the 
address  for  a  location.  Coupled  with  DEPOSIT  and  EXAMINE  keys,  it 
is  thus  possible  to  examine  and  change  any  location  in  any  memory 
while  the  computer  is  halted.  This  process  will,  however,  destroy 
the  previous  contents  of  register  2  of  the  X-memory  (Black  Box 
DMA  address)  ,  since  this  location  is  used  to  specify  the  black  box 
address  of  all  I/O  transfers. 

In  addition  to  the  DEPOSIT  and  EXAMINE  keys,  START,  STOP, 
CONTINUE,  and  SINGLE  STEP  will  be  provided,  as  well  as  EXAMINE 
NEXT  and  DEPOSIT  NEXT.  The  address  for  the  first  instruction  to 
be  executed  following  START  is  located  in  register  2  of  the  X-memory. 

Additional  display  lights  will  include  the  Instruction 
Register  (IR) ,  Input  Buffer  Register  (IBR) ,  and  Output  Buffer 
Register  (OBR) . 

9.  Miscellaneous 

a)  Timing:  It  is  anticipated  that  all  instructions 
except  multiply  will  tako  200  nsec,  and  multiply  will 
require  300  nsec.  These  are  conservative  figures  and, 
as  more  packages  become  available  in  Schottky  TTL,  will 
probably  be  revised  downward. 

b)  Size :  We  expect  that  the  entire  processor  will  require 
5  1/4"  rack  height.  About  300-400  IC's  will  be  required 
(the  multiplier  alone  requires  100  74Sl81's).. 
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10.  Programming  Examples:  In  this  section  we  show  several  programs, 
ranging  from  the  straightforward  to  the  spectacularly  obscure,  which 
illustrate  design  features  of  the  machine. 

a.  Complex  multiplication: 

**  j  b. 

C  \  ^ 

(..<•  -bO  '  i 

£>  *T  Uv  V/ 


Initially,  let  X  and  Y  index  registers  contain  N  and  L  respectively, 
and  the  D  index  register  contain  M.  The  result  is  to  go  in  loc¬ 
ations  M  and  M+l.  Therefore,  initially,  we  have; 


(/-\+  V\  ^  UJ-H  W>c 


/  (  C«-)  S 

j  Cw0  ' 

j  ^  a.  — i 

j  (_  be)  ->  V 

y/  ^  (k  ^  H  t  \ 
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* 
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Note:  "*"  means  enable  index  register  for  indicated  field.  "*1" 


1  ORed  with  the  appropriate  index  register. 
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b*  Second  Order  Difference  Equation:  The  second  order  diff¬ 
erence  equation, 

may  be  implemented  by  setting  aside  two  storage  locations:  Y1  for 
yn-l  '  :Uid  Y2  ^or  ^n-2  '  and  Performing  the  following  sequence; 

i  Ga>>.x  -*  T 
V-> 

^  X-  ^  \  O  j  {*\  2  v-fyS^ey 

8/  1 1>  /  (_  y.  -»  G  •  VI  -*  2-) 

Y\  O-i  V2_ 

A/  /  Cz-  *  -*  ) 

Clearly,  a  coefficient  C  other  than  unity  for  xr  can  be  handled 
simply  by  replacing  the  first  MOVE  by  a  MUL  (  MUL  C,  X,  10  ) . 

It  should  be  noted  that  at  the  end  of  the  sequence,  the  Z  register 
(the  addition  entry  port  to  the  multiplier)  has  been  cleared 
since  the  last  multiply  didn't  store  in  it. 

c*  Sum  and  Difference  of  Two  Buffers:  For  convenience,  we 

assume  the  four  buffers  start  at  locations  1,  101,  201,  and  301, 

and  that  the  length  N  of  the  buffers  is  =  100.  Also,  A.  and  B 

i  i 

can  be  in  "overlapping"  X  and  Y  locations,  and  the  Sum  and  Diff¬ 
erence  could  be  both  in  X  or  both  in  Y,  or  both  in  both  X  and  Y. 


which  is  coded  as, 
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e.  Blt-rever3ed  Counter;  Given  a  register,  COUNT,  containing 
a  number  to  be  incremented  in  bit  reversed  order,  and  another  reg¬ 
ister,  WHERE,  containing  the  number  of  the  bit  position  (LSB=0) 
at  which  the  count  should  start: 
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The  following  program  works  because  in  single-source  instructions, 
if  the  Y  index  control  bit  is  1,  the  Y  index  register  will  be  ORed 
with  the  bits  of  the  instruction  which  are  in  the  same  position  as 
the  8-bit  Y  address  of  double-source  instructions: 
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BITC  (test  bit  and  complement  it  )  will  skip  until  the  counting 
operation  is  complete.  Thus  LOOP  is  the  even  location  reached  on 
counting  done,  and  LOOP+1  contains  the  instruction  that  steps  the 
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test  bit  loeation  to  the  right,  and  jumps  to  the  BITC  test  of  the 
new  selected  bit. 


f.  Division:  The  following  mnemonics  are  assumed: 

MOVLO  shift  in  Ones  from  the  right,  and  shift  out  through 
the  link 

MOVLZ  shift  in  Zeroes  from  the  right,  and  shift  out  through 
the  link 

MOVL  rotate  left  through  the  link 
(  By  addition  of  one  instruction  to  the  routine  (preset  the  link) , 
the  requirements  can  be  reduced  to  only  the  availability  of  MOVL.) 

The  following  is  a  routine  for  the  positive  integer  divide  of 
the  double  precision  number  (A , B)  by  C. 
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Abstract — When  digital  signal  processing  operations  are  imple¬ 
mented  on  a  computer  or  with  special-purpose  hardware,  errors  and 
constraints  due  to  finite  word  length  are  unavoidable.  The  ma:a  cate¬ 
gories  of  finite  register  length  effects  are  errors  due  to  A/D  conver¬ 
sion,  errors  due  to  roundoffs  in  the  arithmetic,  constraints  on  signal 
levels  imposed  by  the  need  to  prevent  overflow,  and  quantization  of 
system  coefficients.  The  effects  of  finite  register  length  on  implemen¬ 
tations  of  linear  recursive  difference  equation  digital  filters,  and  the 
fast  Fourier  transform  (FFT),  are  discussed  in  some  detail.  For  these 
algorithms,  the  differing  quantization  effects  of  fixed  point,  floating 
point,  and  block  floating  point  arithmetic  are  examined  and  com¬ 
pared. 

The  paper  is  intended  primarily  as  a  tutorial  review  of  a  subject 
which  has  received  considerable  attention  over  the  past  few  years. 
The  groundwork  is  set  through  a  discussion  of  the  relationship 
between  the  binary  representation  of  numbers  and  truncation  or 
rounding,  and  a  formulation  of  a  statistical  model  for  arithmetic 
roundoff.  The  analyses  presented  here  arc  intended  to  illustrate 
techniques  of  working  with  particular  models.  Results  of  previous 
work  are  discussed  and  summarized  when  appropriate.  Some  ex¬ 
amples  are  presented  to  Indicate  how  the  results  developed  for 
simple  digital  filters  and  the  FFT  can  be  applied  to  the  analysis  of 
more  complicated  systems  which  use  these  algorithms  as  building 
blocks. 

I.  Introduction 

IN  PRACTICE,  digital  signal  processing  requires  t he  rep¬ 
resentation  of  sequence  values  in  a  binary  format  with  a 
finite  register  length.  The  effect  of  the  finite  word-length 
constraint  manifests  itself  in  several  different  ways.  If  a  se¬ 
quence  to  lie  processed  is  derived  by  sampling  an  analog 
waveform,  then  the  finite  word-length  constraint  requires 
that  the  analog-to-digitai  conversion  produce  only  a  finite 
number  of  values.  This  represents  quantization  of  the  input 
waveform.  Even  when  we  start  with  data  representable  with 
a  finite  word  length,  the  result  of  processing  will  naturally 
lead  to  values  requiring  additional  bits  for  their  representa¬ 
tion.  For  example,  a  6-hit  data  sample  multiplied  by  a  6-hit 
coefficient  results  >n  a  product  which  is  26  hits  long.  If  in  a  re¬ 
cursive  digital  filter  we  do  not  quantize  the  result  of  arith¬ 
metic  operations,  t he  number  of  hits  required  will  increase 
indefinitely,  since  after  the  first  iteration  26  bits  are  required, 
after  the  second  iteration  36  bits  are  required,  etc.  The  effect 
of  quantization  in  sucli  a  context  depends  on  such  factors  as 
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whether  we  are  considering  fixed-point  or  floating  point  arith¬ 
metic,  and  whether  for  fixed-point  arithmetic  wt  are  using  a 
representation  of  numbers  in  terms  of  fractions  or  integers, 
or  perhaps  a  mixture.  We  will  he  treating  the  case  of  fixed- 
point  arithmetic  and  floatiugqioint  arithmetic  separately. 
For  fixed-point  arithmetic,  it  is  natural  in  a  signal  processing 
context  to  consider  a  register  as  representing  a  fixed-point 
fra.tion.  In  this  way  the  product  of  two  nnmlnTs  remains  a 
fraction  and  the  limited  register  length  can  lie  maintained  by 
truncating  or  rounding  the  least  significant  hits.  With  this 
type  of  representation  the  result  of  addition  on  fixed-|>oint 
fractions  need  not  lie  truncated  or  rounded  hut  it  can  increase 
in  magnitude  so  that  the  sum  eventually  is  not  a  fraction. 
This  effect  is  commrnly  referred  to  as  overflow,  and  can  Ik* 
handled  by  requiring  that  the  input  data  he  siifficicntiv  small 
so  that  the  possibility  of  overflow  is  avoided.  In  considering 
floating-point  arithmetic,  dynamic  range  considerations  gen¬ 
erally  can  he  neglected  due  to  the  large  range  of  representable 
numbers,  luit  quantization  is  introduced  both  for  multiplica¬ 
tion  and  for  addition. 

A  third  effect  of  finite  word  length  is  inaccuracies  in  pa¬ 
rameter  values.  While  generally  signal  processing  parameters 
are  initially  specified  with  unlimited  accuracy,  they  can  only 
lie  utilized  with  finite  word  length  This  effect  is  similar  to  the 
effect  which  arises  in  implementing  analog  pr>  cessing  using 
inaccurate  circuit  elements.  There  are  two  possible  approaches 
to  handling  tbe  inaccuracies  in  parameter  values.  One  |>ossi- 
bility  is  to  develop  design  procedures  which  inherently  are 
insensitive  to  parameter  inaccuracies.  An  alternate  is  to  choose 
specifications  which  are  consistent  with  the  limited  register 
length.  There  is  a  certain  amount  that  is  understood  about 
the  effect  of  inaccuracies  in  parameter  values,  hut  for  the 
most  part  present  results  lead  to  guidelines  rather  than  hard 
design  or  analytical  strategies. 

In  the  following  discussion  the  relationship  lietween  the 
binary  representation  of  numbers  and  truncation  or  rounding 
is  discussed  and  a  statistical  model  for  arithmetic  roundoff  is 
presented.  This  statistical  model  is  then  applied  to  the  anal¬ 
ysis  of  fixcd-|>oint  and  floatingqioiiit  rounding  errors  in 
digital  filters.  The  analysis  includes  a  consideration  of  the 
effect  of  dynamic  range  in  developing  and  comparing  signal 
to-noise  ratios  for  fixed-point  and  floating-point  filters.  It 
is  not  always  possible  to  treat  the  effects  of  arithmetic  round 
off  in  terms  of  a  simple  statistical  model  Some  approaches 
and  results  are  available  in  the  literature  on  the  limit  cycle 
Ircbavior  of  digital  filters  due  to  arithmetic  roundoff,  anil  a 
discussion  of  some  of  these  results  is  included. 

For  the  analysis  of  arithmetic  roundoff  in  computation  of 
the  discrete  Fourier  transform  using  the  fast  Fourier  trans- 
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form  (FFT)  algorithm  a  statistical  model  is  used.  With  this 
model  the  signal-to-noise  ratio  is  developed  and  compared 
for  fixed-point  and  floating-point  arithmetic. 

While  for  any  given  filter  configuration  or  spectral  anal¬ 
ysis  problem  it  can  be  difficult  to  carry  out  a  detailed  anal¬ 
ysis  of  the  effects  of  finite  register  length  there  are  a  number 
of  general  guidelines  that  can  lie  distilled  from  the  results  pre¬ 
sented  here.  In  Section  IV  some  examples  and  guidelines  are 
presented  for  filters  implemented  with  fixed-point  arithmetic 
and  with  floating-|ioint  arithmetic  as  well  as  for  filters  im¬ 
plemented  with  the  FFT. 

This  paper  is  intended  primarily  as  a  tutorial  review  of  a 
subject  which  has  received  considerable  attention  over  the 
past  few  years.  The  analyses  which  are  presented  here  are 
selected  to  illustrate  techniques  of  working  with  particular 
models.  Previous  work  is  freely  referenced,  discussed,  and 
borrowed  from. 

II.  Number  Representation  and  Its 
Effect  on  Quantization 

A.  Fixed-Point  and  Floating-Point  Numbers 

The  manner  in  which  finite  word-length  effects  are  mani¬ 
fested  is  closely  tied  to  the  way  in  which  numbers  are  repre¬ 
sented. 

Digital  computers  and  special  purpose  digital  hardware 
for  the  most  part  use  a  number  representation  with  a  radix  of 
2,  i.e.,  a  binary  representation.  Therefore,  a  number  is  repre¬ 
sented  by  a  .quence  of  binary  digits  which  are  either  zero 
or  unity.  Just  a  decimal  number  is  represented  as  a  string 
of  decimal  dig.ta  with  a  decimal  point  dividing  the  integer 
part  from  the  fractional  part,  the  sequence  of  binary  digits 
is  divided  by  a  binary  |w>int  into  those  representing  the  in¬ 
teger  part  of  the  number  and  those  representing  the  frac¬ 
tional  part.  Thus  if  A  denotes  'he  location  of  the  binary  |>oint, 
the  binary  number  IOOI4OUO  has  the  decimal  value  of  (I  X2* 
+0  X  2s + 0  X  2 1 +  1  X  2") +  (0  X  2-' + 1  X  2-*+ 1  X  2-» +0  X  2^ '). 
This  representation  always  corresponds  to  a  positive  number. 

The  manner  in  which  arithmetic  is  implemented  in  a 
digital  computer  or  in  a  special  purpose  hardware  de|iends  on 
where  in  the  register  the  binary  point  is  located.  For  fixed- 
point  arithmetic,  the  implementation  is  based  on  the  assump¬ 
tion  that  the  location  of  the  binary  [mint  is  fixed.  The  manner 
in  which  addition  is  carried  out  will  not  depend  on  the  location 
of  the  binary  poi nt  for  fixed-point  arithmetic  as  long  as  it  is 
thi  same  for  every  register.  For  multiplication,  however,  the 
location  of  the  binary  |«>int  must  be  known.  For  example, 
consider  the  product  of  the  two  4-bit  numliers  1001a  and 
0011a.  In  general,  of  course,  the  product  of  two  6-bit  numlters 
will  lie  26  bits  long.  The  8-bit  product  of  the  above  numlier  is 
00011011a-  If,  on  the  other  hand,  we  consider  the  4-bit  frac¬ 
tions  a  1001  and  aOOII,  then  the  8-bit  product  is  a00011011. 
In  digital  filtering  applications,  it  is  usually  necessary  to  ap¬ 
proximate  the  26-bit  product  of  two  6-bit  numbers  by  a  6- 
bit  result.  In  integer  arithmetic  this  is  difficult.  With  frac¬ 
tional  arithmetic,  on  the  other  hand,  this  can  lie  accomplished 
by  truncating  or  rounding  to  the  most  significant  6  bits.  For 
multiplication  with  fractions,  overflow  can  never  occur  since 
the  product  of  two  fractions  is  a  fraction.  Thus  for  the  4-bit 
example  previously  mentioned,  the  product  aOOOUOH  ran  lie 
approximated  by  aOOOI  (truncation)  or  a0010  (rounding). 

An  alternative  to  fixedqioint  arithmetic  is  a  floating-point 
representation.  In  this  case,  a  positive  numlier  Fisrenresented 
as  F=P.\I,  where  .1/,  the  mantissa,  is  a  fraction  lietween  1  2 
and  1,  and  c,  the  characteristic,  can  lie  either  positive  or 


negative.  The  product  of  two  floating-point  numbers  is 
carried  out  by  multiplying  the  mantissa  as  fixed-point  frac¬ 
tions  and  adding  the  characteristics.  Since  the  prtiduct  of  the 
mantissas  will  be  lietween  1/4  and  1,  a  normalization  of  the 
mantissa  and  corresponding  adjustment  of  the  characteristic 
may  be  necessary.  The  sum  of  two  floating-point  numbers  is 
carried  out  by  scaling  the  mantissas  of  the  smaller  numlier  to 
the  right  until  the  characteristics  of  the  two  numbers  are 
equal  and  then  adding  the  mantissas.  For  example,  consider 
the  sum  of  Ft  and  Ft  with  Fi  =  4  and  F;  =  5/4.  Then  in  floating¬ 
point  notation,  F,  =  2**3/i,  and  Ft  =  P'M,  with 

C\  =  11a  (  =  3  decimal) 

M\  =  al000  (  =  0.5  decimal) 

«j  =  1a  (=1  decimal) 

Mi  —  al010  (  =  5/8  decimal). 

In  order  to  carry  out  the  addition,  ct  must  be  changed  to  equal 
Ci  and  Mi  must  be  adjusted  accordingly.  Thus  first  the  repre¬ 
sentation  of  Ft  is  changed  to  Fj  =  2f,i/s  with 

it  =  Ha 

=  a00101 

in  which  case  the  mantissas  can  now  be  added.  The  resulting 
sum  is  F=  P M  with  c  =  11a  and  .1/  =  a  101 01.  In  this  case  the 
sum  of  Mi  and  ifi  is  a  fraction  lietween  1  /2  and  1  and  therefore 
no  further  adjustment  of  c  has  to  lie  carried  out.  In  a  more 
general  case,  the  sum  may  not  lie  in  that  range,  and  conse¬ 
quently,  c  would  lie  adjusted  to  bring  the  mantissa  into  the 
pro|ier  range.  From  this  example  it  should  lie  clear  that  in 
general  with  floating-point  arithmetic,  the  mantissa  can  ex¬ 
ceed  the  register  length  and  must  therefore  lie  truncated  or 
rounded  for  both  addition  and  multiplication  whereas  this  is 
only  necessary  for  multiplication  in  the  fixed-point  case.  On 
the  other  hand,  if  the  result  of  addition  in  the  fixed-point 
rase  exceeds  the  register  length,  truncation  or  rounding 
will  not  help,  i.e..  the  dynamic  range  has  lieen  exceeded.  Thus 
while  floating  point  introduces  error  due  to  arithmetic  round¬ 
off,  it  provides  much  greati  r  dynamic  range  than  fixed  point. 
As  we  will  see  later,  both  f  these  effects  must  lie  considered 
when  comparing  fixed-point  and  floating-point  realizations 
of  digital  filters. 

B.  Representation  of  Negative  Numbers 

There  are  three  common  means  used  for  representing 
fixed-point  negative  numliers.  The  first,  and  most  fa>..i!i.:. ,  is 
sign  and  magnitude,  i.e.,  the  magnitude  (which  is  of  course 
positive)  is  represented  as  a  binary  number  and  the  sign  is 
represented  by  the  leading  binary  digit  which,  if  0  corresponds 
to  a  +  and  if  1  corres[Hinds  to  a  —  (or  vice  versa).  Thus  for 
example,  in  sign  and  magnitude  0a0011  represents  3  16  and 
1  a tX)  1 1  represents  —3  16.  Two  other  related  representations 
of  negative  numliers  are  often  referred  to  as  one's-comple- 
ment  and  two's-complement  representations.  Considering  all 
numliers  to  lie  fractions,  a  positive  numlier  is  represented  as 
liefore.  For  two's  complement  representation  negative  num¬ 
lier  is  represented  by  2.0  minus  its  magnitude.  For  example 
—  (OaOUO)  in  sign  and  magnitude  is  represented  as  l^lOtO 
in  two's-complement  since  IOaOOO  —  Oa0100  =  IaIOIO.  For 
one's-complement,  the  negative  num'ier  is  represented  by 
subtracting  the  magnitude  from  the  largest  numlier  repre¬ 
sentable  in  the  register.  Thus  —(OaOUO)  is  represented  by 
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(latlll)  —  (OiOllO)  =  ljlOOl.  One’s  complement  representa¬ 
tion  is  equivalent  to  representing  a  negative  numlier  by  the 
bit-by-bit  complement  of  its  magnitude.  The  choice  of  repre¬ 
sentation  for  negative  numbers  in  a  particular  system  is 
usually  based  almost  entirely  on  hardware  considerations. 

For  the  representation  of  negative  floating- point  numbers 
there  are  a  variety  of  conventions  that  have  been  used.  In 
this  paper  we  will  consider  the  sign  of  the  numlier  to  be  asso¬ 
ciated  with  the  mantissa  so  that  the  mantissa  is  a  signed 
fraction.  The  representation  of  this  signed  fraction  can  of 
course  be  in  sign  and  magnitude,  one’s-complemcnt  of  two's- 
complement  notation. 

C.  A  Model  for  Arithmetic  Roundoff 

In  formulating  a  model  for  arithmetic  roundoff,  we  shall 
consider  both  fixed-point  numbers  and  mantissas  of  floating¬ 
point  numbers  to  be  represented  as  64- 1-bit  binary  fractions, 
with  the  binary  point  just  to  the  right  of  the  highest  order  bit 
(or  sign  bit).  This  convention  represents  no  loss  of  generality, 
and  its  convenience  has  been  alluded  to  above.  The  numerical 
value  (for  positive  numbers)  of  a  one  in  the  least  significant 
bit  is  2~*,  and  this  quantity  can  lie  referred  to  as  the  width  of 
quantization. 

As  indicated  previously,  the  effect  of  finite  register  length 
on  the  result  of  arithmetic  operations  de|iends  on  whether 
fixed-point  or  floating-point  arithmetic  is  used,  and  how  nega¬ 
tive  numbers  are  represented.  Let  us  consider  first  the  effect 
of  truncation  and  rounuing  in  the  Axed-fmint  case.  For  sign 
and  magnitude,  one's-complement  and  two’s-complement, 
the  representation  of  |Misitive  numliers  is  identical  and,  con¬ 
sequently,  so  is  the  effect  of  truncation  and  rounding  If  Et 
denotes  the  error  due  to  truncation,  i.e.,  the  value  after  trun¬ 
cation  minus  the  value  liefore  truncation,  this  error  will  al¬ 
ways  lie  negative  for  positive  numliers.  That  is,  the  effect  of 
truncation  is  to  reduce  the  value  of  the  numliers.  More 
specifically,  if  6,  denotes  the  numtier  of  bits  (exclusive  of  sign) 
after  truncation,  and  6|  denotes  the  numlier  of  bits  I  efore 
truncation,  then  the  result  satisfies  0>£j->  — (2“*’— 2“*  ). 

With  sign  and  magnitude  representation  of  negative  n  mi¬ 
llers,  truncation  reduces  the  magnitude  of  the  number  and  tto 
error  Et  satisfies  0<£r<(2-**—  2-*').  For  a  two's-comple- 
ment  negative  number  represented  by  the  bit  string  la,  <i,. 
n2,  •  •  •  ,  as,,  the  magnitude  is  given  by 

Mi  =  2.0  -  r, 


where 


X\  =  1  4-  52 

i-i 

Truncation  to  6-.  bits  (6j<6i)  produces  the  bit  string  la,  at, 
Uj,  •  •  •  .  at,,  where  now  the  magnitude  is 


M,  =  2.0  -  xt 


with 


*« 

**  =  1  +  £  th2-*. 

>-l 


The  change  in  magnitude  is 


ooo 


0(a)  Oil) 


Rounding 

-•(•***«  «(•>-» 


Trgncotion 

(2*«  complainant) 

0  2  0(a)-a  >2-* 


Truncation 

(<’•  complainant  and  sign 
and  mognituda 
0Z0(«)-i»2-“  ;  i>0 
0S0(E]-.<Z-»,  i«0 


Fig.  I.  Transfer  characteristics  for  rounding  and  truncation. 


and  it  is  easily  seen  that 

0  <  AM  <  2-‘»  -  2-»*. 

Hence  the  effect  of  truncation  for  two’s-complement  negative 
numbers  is  to  increase  the  magnitude  of  the  negative  number: 
the  truncation  error  is  negative,  and  satisfies  0  >£r>  —  (2“** 
-2-»‘). 

For  a  one’s-complement  negative  number  represented  by 
the  bit  string  1it  at,  o2,  •  •  •  ,  a*,,  the  magnitude  is  given  by 

Mi  =  2.0  -  2-*'  -  Xi 

and  truncation  to  6;  bits  yields  a  magnitude 
Ms  =  2.0  -  2-**  - 

where  .Vi  and  at  are  as  defined  above.  The  change  in  magnitude 

is 


b  I 

AM  =  Mi  -  M I  =  £  <1.2-'  -  (2-*  -  2-‘<) 

and  now 

0  >  AM  >  -  (2-»-  -  2-*'). 

Hence  the  effect  of  truncation  for  one’s  complement  negative 
numliers  is  to  decrease  the  magnitude  of  the  negative  num 
tier:  the  truncation  error  is  positive,  and  satisfies  0<£r<2"®* 
-  2*'. 

The  effect  of  rounding,  of  course,  will  be  the  same  inde¬ 
pendent  of  how  negative  numliers  a-  represented  and  the 
rounding  error  will  always  be  greater  than  or  equal  to 
(  —  1/2)2“*  and  less  than  or  equal  to  (•+•  l/2)2~*.  The  effect  of 
truncation  and  rounding  for  the  fixed-|Kiint  case  is  immarized 
in  Fig.  1,  where  at  represents  the  value  before  truncation  or 
rounding  and  Q(x)  represents  the  value  after.  In  the  figure  it 
is  assumed  that  x  can  take  on  a  cont'nuous  range  of  values, 
corresponding  to  6<=  «  in  the  discussii  n  above,  and  that  the 
quantized  word  length  is  6  bits  plus  sign. 

For  the  case  of  floating-point  arithmetic,  the  effect  of 
truncation  or  rounding  is  reflected  only  in  the  mantissa.  It  is 
convenient  in  the  floating-point  case  to  descrilie  the  error  in  a 
multiplicative  sense  rather  than  in  an  additive  sense  as  is 
done  in  fixed-point  arithmetic.  In  other  words,  for  a  floating¬ 
point  word,  if  x  represents  the  value  liefore  truncation  or 
rounding  and  Q(x)  represents  the  value  after,  then  we  express 
Q(x)  as  equal  to  art  1  4-*)  For  the  case  of  rounding,  for  example, 
the  error  in  the  mantissa  is  between  ±2~*  2.  and  conse¬ 
quently  the  error  in  the  value  of  the  floating-point  word  is 


AM  =  Mi  -  Mi 


Z 


1-kjF  1 


<.,2-' 


2“*  2“* 

-2* - <  Q(x)  —  x  <  2* - 

2  2 


SO 


«<  I 
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PIE,) 


(a) 


p<Et> 

|  2  s  Co<n|Samsnt  Truncation 


-V 


(b) 


Fig  2.  (a)  Probability  density  function  for  rounding  noise,  (b)  Proba¬ 

bility  density  function  for  noise  due  to  two's-complement  truncation. 


or,  since Q(x)—x  =  tx 

2-»  2-* 

—  2* - <  s.r  <  2 - 

2  2 

and  since  2e_l  <x<2',  we  can  write  that  (or  the  case  of  round¬ 
ing  —  2-*<s<2-t.  In  a  similar  manner  we  can  show  that  for 
one’s-complement  and  for  sign  and  magnitude  truncation 
0>s>  —  2-2~*.  For  two’s-complement  truncation 

<)>.>-  2-2-*,  x  >  0 

0<«<2-2-*,  *<0. 

D.  Statistica I  Model  of  Arithmetic  Roundoff 

A  convenient  means  for  analyzing  the  effect  of  quantiza¬ 
tion  :s  to  represent  the  error  statistically  ( 1  ].  [2].  In  par¬ 
ticular.  for  the  case  of  fixed-point  arithmetic  and  rounding  h  r 
is  represented  as  a  random  variable  with  a  probability  density 
shown  in  Fig.  2(a).  For  the  case  of  two’s-complement  trunca¬ 
tion,  the  probability  density  is  shown  in  Fig.  2(b). 

In  each  of  these  cases,  the  assumption  is  that  the  random 
variable  Et  is  independent  of  x.  For  one’s  complement  and 
sign  magnitude  truncation,  this  assumption  cannot  lie  made 
since  the  mean  value  of  the  error  is  directly  correlated  with 
the  sign  of  .v.  In  the  analysis  that  follows  for  fixed-point  arith 
metic,  the  discussion  is  phrased  in  terms  of  rounding.  The  re¬ 
sults  are  easily  modified  for  two's-complement  trum  .  lion.  In 
|>artirular  the  variance  of  the  noise  is  identical  for  both  cases. 
However,  for  rounding  the  noise  is  zero  mean  and  for  two’s- 
complement  truncation  it  is  not  zero  mean. 

For  the  ffoating-point  case,  the  parameter!  is  con-idered 
to  be  a  random  variable  which  is  independent  of  x.  In  that 
case  the  assumption  of  independence  is  reasonable  for  round¬ 
ing,  sign  and  magnitude  truncation,  and  one’s-complement 
truncation,  out  not  for  two's-complement  truncation.  The 
random  variables  is  bounded  by  —  2~*<s<2_*.  We  will  gen¬ 
erally  assumes  to  tie  uniformly  distributed  in  this  range  with  a 
variance  <r,*»  (1  3)2~*.  Empirical  work  has  shown  that  the 
distribution  is  not  quite  uniform  so  that  while  »,*  is  promo¬ 
tional  to  2  *.  the  constant  of  proportionality  is  siigitii,  '»ss 
than  l/i.  However,  the  interpretation  of  the  results  depends 
primarily  on  the  proportionality  to  2  * 

III.  Finite  Register  Length  Effects 
for  Digital  Filters  |J] 

A.  Introduction 

The  basic  arithmetic  operations  involved  in  implementa¬ 
tion  of  a  digital  filter  are  multiplication  by  a  constant  and 


addition.  For  fixed-point  arithmetic,  roundoff  is  introduced 
only  after  the  multiplication.  Because  of  the  possibility  of 
overflow  due  to  addition,  there  is  a  dynamic  range  limitation 
in  fixed-point  fillers.  In  contrast,  floating-point  filter  imple¬ 
mentation  has  a  much  less  severe  dynamic  range  constraint, 
although  arithmetic  roundoff  is  introduced  due  to  both  multi¬ 
plication  and  addition.  In  the  next  sections  we  will  first  de¬ 
velop  the  statistical  analysis  of  arithmetic  roundoff  for  fixed- 
point  fillers  including  dynamic  range  considerations.  This  is 
followed  by  a  statistical  analysis  for  floating-point  arithmetic 
and  a  discussion  of  zero  input  limit  cycle  behavior  for  fixed 
point  arithmetic. 

B.  Statistical  Analysis  of  Fixed-Point  Errors  in  a  Digital 
Filler  (4).  |i] 

In  many  situations  it  is  reasonable  to  model  the  effect  of 
rounding  in  a  digital  filter  by  a  simple  statistical  nu>del.  The 
approach  is  to  model  the  effect  of  the  rounding  at  each  multi¬ 
plier  by  a  white-noise  source  uniformly  distributed  in  anifsli- 
tude  between  plus  and  minus  (1  2)2  6  Each  of  the  noise 
sources  is  assumed  to  be  linearly  inde|>endent  of  each  other 
and  of  the  input.  Experimentally  these  assumptions  have 
been  justified  for  a  broad  class  of  inputs  including  random 
signals,  speech,  etc.  The  model  is  clearly  not  valid  for  certain 
inputs,  such  as  constant  inputs.  If  the  impulse  response  from 
the  Fth  noise  source  to  the  output  is  **(«)  then  the  steady- 
state  output  noise  variance  due  to  the  k th  noise  source  is 


*«*’“».*  E  **’( n)  (1) 

where  12)2  *  Since  all  the  noise  sources  are  as¬ 

sumed  to  be  uncorrelated,  the  total  output  noise  is 

*.*  -  E  *-’•  (2) 

s 


For  example,  if  we  consider  the  first-order  filter  in  Fig.  I  one 
noise  source  is  introduced.  In  this  case,  the  impulse  response 
from  the  noise  source  input  to  the  output  is  *(«)=«"«  ,(*) 
where  denotes  a  unit  step  sequence,  so  that 


*» 


> 


1 

12 


0) 


F’or  a  secoi.,1  order  filter  with  o  te  i  omplex  pole  pair  there  are 
two  noise  sources  as  indicated  ii,  Fig  4  The  resulting  output 
noise  is 


2  / 1  +  r*  I  \ 

f.‘  -  —  2 - 1. 

12  \1  -  r*  r*+  I  -  2r’  cos  29  / 


(4) 


C.  Dynamic  Range  Considerations  for  Fixed-Point  Filters 

As  indicated  previously,  the  possibility  of  overflow  must 
be  considered  in  the  implementation  of  digital  filters  with 
fixedqsiint  arithmeti  .  With  the  convention  that  each  fixed- 
point  register  represents  a  signed  fraction,  each  node  in  the 
filter  must  be  constrained  to  maintain  a  magnitude  less  than 
unity  in  order  to  avoid  overflow  l-rtting  *(•)  denote  the  filter 
input  and  y»(n)  and  kt(n)  denote  the  output  and  unit  sample 
response  few  the  kth  node  in  the  filter,  then 


?*(«)  “  E  *»(*)*(■  “  •’)■  (5) 

Ml 

If  r..,  denotes  the  maximum  of  the  absolute  value  of  the  in- 
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that  the  dynamic  range  of  the  registers  is  not  exceeded.  If  we 
consider  the  input  sequence  to  be  uniformly  distributed  white 
noise,  the  resulting  output  noise-to-signal  ratio  will  be 


Fig.  3.  Noisy  first-order  filter  (fixed  point). 


While  it  is  difficult  to  evaluate  this  expression  exactly,  it  is 
|M>ssible  to  obtain  an  upper  and  lower  bound.  Since  //n|  is 

the  largest  |>ossible  output  obtainable  with  an  input  that 
never  exceeds  unity,  it  must  be  larger  than  the  response  of 
the  second-order  filter  to  a  sinusoid  of  unity  amplitude  at  the 
resonant  frequency.  With  this  consideration,  we  can  write  that 


since  the  right-hand  side  of  this  inequality  is  the  gain  at  reso¬ 
nance.  Furthermore, 


Fig.  4.  Noisy  second-order  filter  (fixed  point). 


Thus,  since  we  require  that  |  ,v*(i»)|  <1,  (6)  requires  that 


Equation  (7)  thus  provides  an  upper  bound  on  the  maximum 
value  of  the  input  to  insure  that  no  overflow  occurs  in  the  Fth 
n<»de.  For  a  general  input  (7)  in  fact  provides  a  least  upper 
bound,  i.e.,  if  the  maximum  value  of  the  input  exceeds  the 
bound,  overflow  can  occur.  This  is  a  consequence  of  the  fact 
that  equality  can  l>e  achieved  in  (6)  with  a  sequence  x(if)  for 
which  at  n  =  not  x(n0-r)=  [sgn  A*(r)]  for  r  =  0  to  »  (W  here 
sgn  (x) «  1  (or  x>0  and  sgn  (x)  =  —1  for  x<0.)  Thus  *  the 
most  general  case,  (7)  is  required  to  guarantee  that  no  er- 
flow  occurs.  The  condition  in  (7)  would  generally  l*»  sat  sfied 
by  applying  attenuation  to  the  signal  at  the  filter  input. 

If  we  assume,  for  example,  that  the  input  x(n)  is  a  white- 
noise  sequence  with  a  uniform  amplitude  distribution,  we 
would  choose  for  the  rase  of  the  first-order  filter  a  maximal 
input  amplitude  of  fl— a).  For  this  case,  if  a,1  denotes  the 
variance  of  the  input  signal,  and  ej*  denotes  the  variance  of 
the  output  signal,  then 


For  both  the  first-  and  second-order  filter  an  expression  for 
the  noise-to-signal  ratio  can  lie  obtained  which  provides  some 
insight  into  the  behavior  of  the  noi**-to-signal  ratio  as  the 
poles  approach  the  unit  circle.  For  the  first-order  filter  let 
4  =  I— a  so  that  as  4  — *0.  the  pole  approaches  the  unit  circle. 
Then  in  terms  of  6,  the  noise-to-signal  ratio  foi  thy  first-order 
filter  is 


For  the  second-order  filter,  let  4*  1—  r  so  that,  again,  as  6—*) 
the  |K»les  approach  the  unit  circle.  Then  if  we  assume  Cat 
4«1.  we  can  approximate  (1  +r*—  2r  cos  20)  as 

(1  +  r*  -  It  cos  20)  i  4  sin2  0  +  a2  (15) 

which  for  4  sin*0  large  compared  with  4*  we  will  approximate 
as  4  sin*  0.  Consequently,  incorporating  this  approximation. 


-nr  sin- ir  i  o-sin-s 

Thus  we  observe  that  the  n<»iae-to>signal  ratio  as  considered 
thus  far  can  l»e  considered  to  1**  proportional  to  2  *  4s.  We 
note  fr»m  this  dependent  that  if  4  is  halved,  then  to  maintain 
the  same  noise-to-signal  ratio  b  must  be  increased  by  1,  i.e., 
one  Im t  must  l»e  added  to  the  register  length  Tliia dependence 
provides  a  convenient  basis  for  com}narison  of  different  over¬ 
flow  strategies  and  different  kinds  of  arithmetic. 

In  a  similar  manner  we  can  derive  a  noise-to-signal  ratio  for  In  the  almve  analysis,  the  filter  input  was  assumed  to  be 

the  second-order  filter  shown  in  Fig  4.  As  in  the  first-order  uniformly  distributed  white  noise.  As  4  approaches  rero  the 
case,  we  restrict  the  maximum  input  in  order  to  guarantee  frequency  response  of  both  the  first-  and  second-order  filter 


F.»r  this  example,  we  can  then  compute  a  noise -to- s tgnai  ratio 
as  the  ratio  o,}  with  the  result 
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becomes  more  selective  so  that  more  and  more  of  the  input 
energy  is  out  of  band.  An  alternative  basis  for  determining 
the  noise-to-signal  ratio  is  for  an  input  which  is  sinusoidal. 
For  this  choice  of  inputs,  of  course,  we  would  not  use  the 
general  condition  of  (7)  to  avoid  overflow  since  we  can  deter¬ 
mine  exactly  the  maximum  allowable  input  amplitude  as  a 
function  of  the  filter  parameters. 

In  particular,  if  the  input  is  of  the  form  x(n)  cos  n* 

then  the  steady-state  output  is  of  the  form  y(n)  =>w,  cos 
(nf+ifi).  To  prevent  overflow,  yMI  must  be  less  than  unity 
and  to  maximize  the  output  signal  energy,  y„,„  is  chosen  to  be 
as  large  as  possible.  Thus  the  maximum  noise-to-signal  ratio 
is  obtained  when  v„,„  is  chosen  so  that  y(n)  =cos 
Note  that  in  order  to  choose  x™,  in  this  way,  the  frequency  of 
the  input  signal  must  be  known.  For  an  input  sinusoid  of  un¬ 
known  frequency  rm„  must  tie  attenuated  so  that  overflow 
will  not  occur  even  in  the  worst  case,  where  the  frequency  of 
the  input  coincides  with  the  peak  gain  in  the  filter's  transfer 
function. 

For  fixed-point  filters,  within  the  validity  of  the  statistical 
model  for  roundoff  error,  the  output  noise  is  independent  of 
the  form  and  amplitude  of  the  input  signal.  Thus  for  this 
choice  of  inputs,  the  noise-to-signal  ratio  obtained  for  the 
first-order  filter  is 


O, 


> 


(17) 


is  extremely  small.  Furthermore,  for  many  filters  it  is  difficult 
to  compute  the  sum  in  (7).  Jackson  [7]  has  formulated  the 
dynamic  range  constrai  ns  on  fixed-point  digital  filters  in 
terms  of  L,  norms.  In  particular,  let  Y(u),  AT (to),  and  H(u) 
denote  the  Fourier  transforms  of  the  filter  output,  input,  and 
system  impulse  response,  respectively.  Then  it  can  be  shown 
in  general  that 

I  ?(«)  I  sil*IM|x||ti/#  +  i/,-i  (2i) 

where  ||//||,  and  ||Jf||,  are  the  Lp  norm  and  L,  norm  of  H(oi) 
and  -Y(«),  respectively,  where  these  norms  are  defined  as 

Ml,= 

and 

-J  |  *(«)!•  do, J  . 

hor  example,  with  H(w)  chosen  as  unity,  i  consequence  of 
(21)  is  that 

l*(«)|  <||-V||f,  all  9  >  1. 

As  another  --onsequence,  if  we  choose  p=  1,  9  =  oo,  and  use 
the  fact  tha  the  /-«,  norm  of  j  .Y (uy  j  is  the  maximum  value 
of  |  -Y (tt>)|  then  we  obtain  the  statement  that 


If.  as  before,  we  let  a  =  1  —  4,  then  for  4«1 


1  2~u 
<r„2  ~  48  4 


(18) 


Thus  in  this  case,  the  noise-to-signal  ratio  is  proportional  to 
i/S  rather  than  1  /42  so  'hat  if  4  is  multiplied  by  1/4  and  the 
register  length  is  increased  by  one  bit,  the  noise-to-signal 
ratio  will  remain  constant.  We  can  consider  the  second- 
order  case  in  a  similar  manner.  Again  for  a  sinusoidal  input, 
the  output  with  maximum  amplitude  has  the  form  y(n)  = 
cos  (nj+jr)  so  that  the  noise-to-signal  ratio  in  this  case  is 

=  1  2-n/i±!!  1 

<r,*  12  \1  -  r*  I  +  r*  -  2r2cos  2S 

Again,  choosing  r=  1  — 4  with  4«1, 


<r»* 

aJ 


2-«* 

44  sin2  9 


(20) 


so  that,  as  with  the  first-order  filter,  the  noise-to-signal  ratio 
is  proportional  to  1/4  rather  than  1  4*.  The  comparison  in  the 
noise-to-signal  ratio  for  a  white-noise  input  and  a  sinusoidal 
input  serves  to  illustrate  the  dependence  of  the  effect  of  dy¬ 
namic  range  considerations  on  the  particular  form  of  the  in¬ 
put.  In  some  sense,  the  two  cases  considered  represent  ex¬ 
tremes.  As  the  input  becomes  more  confined  to  a  known  nar¬ 
row  band  of  frequencies  the  alx.ve  analysis  with  a  sinusoidal 
input  would  lie  more  representative,  and  as  the  input  becomes 
more  wide  band  the  above  analysis  with  a  white-noise  input 
is  more  representative. 

In  the  above  discussion,  the  noise-to-signal  ratio  for  the 
case  of  white-noise  input  was  derived  on  thi  basis  that  over- 


|  y(n)  I  <  max  [  |  .Y (w)  j  j  —  f  |  //(«)  |  dw. 

2w  J 

As  an  alternative,  with  p  =  2,  9  =  2, 

[1  r r  "l1/* 

-J  \hm\uu\ 

To  prevent  overflow  in  the  output  we  require  that 
|  y(«)|  <1  and  to  insure  this  from  (21)  we  will  require  that 
!!7/|UI-Y||,<l.  Consequently,  the  input  must  be  scaled  in 
such  a  way  that 

11*11.  <  I/Mf-  (22) 

This  condition  is  somewhat  less  general  than  (7)  but  in  many 
cases  is  easier  to  apply.  According  to  (22)  with  p  =  2,  ?  =  2, 
the  condition  is  in  terms  of  the  energy  in  the  input  signal  and 
the  energy  in  the  system  impulse  response.  For  9=1,  p  =  «, 
(22)  provides  a  bound  in  terms  of  the  (leak  value  of  the  mag¬ 
nitude  of  the  transfer  function,  which  is  perhaps  most  ap¬ 
propriate  for  a  sinusoidal  input. 

For  the  case  of  a  random  input  (21)  cannot  he  applied 
since  the  input  and  output  do  not  have  Fourier  transforms. 
In  this  case  the  corresponding  condition  is  phrased  in  terms  of 
*»(><)  the  autocorrelation  function  of  the  output,  <!>„(< a)  the 
power  density  spectrum  of  the  input,  and  //(&>)  the  magnitude 
of  the  system  function.  In  particular,  the  inequality  corre¬ 
sponding  to  (21)  is 

♦»(«)  <  MUM.  (23a) 

or  equivalently 


flow  must  lie  avoided.  In  a  practical  case,  a  scaling  of  the  in-  t  \  *■  \\  //!  2  IU 

put  on  the  basis  of  (7)  can  be  considered  to  be  somewhat  pessi-  (2.1b) 

mistic  since  the  probability  of  equality  being  attained  in  (7)  Since,  if  the  input  is  zero  mean.  *„„(())  =»,*  it  follows  that 


f 
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<  MlllMr 

(24) 

Two  particular 

cases  of  interest  are  p=l,  q=  » 

and  «, 

9”  1  so  that 

(25) 

and 

<  ||*|i:i|4>„!|,. 

(26) 

As  Jackson  points  out,  (25)  implies  the  most  stringent  condi¬ 
tion  on  the  input  spectrum  ♦«(«)  whereas  (26)  implies  the 
most  stringent  condition  on  the  transfer  function.  From  (25), 
if  the  input  spectrum  is  white  so  that  ♦,.,(«)  =ox*  for  allot,  then 

a,  <  a,  ||//||.  (27) 

with  the  input  sequence  Gaussian,  then,  the  output  will  over¬ 
flow  no  more  often  than  the  input  overflows  if 

M*  <  1.  (28) 

More  generally,  (27)  provides  a  basis  for  choosing  the  input 
variance  to  control  the  maximum  percentage  of  time  that  the 
output  can  overflow. 

D.  Statistical  Analysis  of  Roundoff  Errors  with  Floating-Point 
Arithmetic 

For  the  case  of  floating-point  arithmetic,  noise  is  introduced 
due  both  to  the  adds  and  the  multiplies.  In  analyzing  the 
effect  of  floating-point  roundoff  the  effect  of  rounding  will  lie 
represented  multiplicatively  so  that  if  [.v]  denotes  rounding  of 
the  mant  ssa  in  a  floating-point  number,  then 

[..]  =  *(1  -I-  «).  (29) 

To  illustrate  the  analysis  of  roundoff  errors  with  floating-point 
■'rithmetic  let  us  consider  a  first-order  filter.  Let  win)  denote 
the  ideal  response  of  the  filter,  that  is,  the  response  with  no 
roundoff  noise  and  let  yin)  denote  the  response  of  the  filter  in 
the  presence  of  roundoff  noise.  Then  following  Liu  and  Kaneko 
[8]  we  can  write  that 

w(n)  =  aw(n  —  1)  +  x(n)  (30) 

yin)  =  |ay(«  -  1)(1  +  *„)  +  .v(»i)](t  4-  £„)•  (31) 

We  assume  that  «„  and  £„  are  uniformly  distributed  between 
—  2~*  and  2“*.  are  uncorrelated  from  iteration  to  iteration,  are 
independent  of  each  other,  and  also  are  independent  of  the 
»  signal.  Letting  E(n)  represent  the  error  in  the  output,  so  that 
Ein)  =  yin)— win),  we  can  write  from  the  aliove  two  equations 
that 

£(»)  -  aF-in  -  1) 

=  awin  -  l)(e.  +  {„)  4-  .r(w)f,  =  u(n)  (32) 

where  we  have  neglected  second-order  terms  in  t.  {,  and  E. 
Since*  and  £  are  statistically  inde|iendent  of  x,  and  of  win  —  1), 
the  term  uin)  is  easily  shown  to  lie  a  white-noise  sequence.  Its 
variance,  of  course,  depends  on  the  excitation  jt(ii).  The 
derivation  of  (32)  with  the  second-order  terms  neglected  cor¬ 
responds  to  representing  the  roundoff  noise  as  an  additive 
noise  source  that  is  statistically  independent  of  the  signal  but 
whose  variance  depends  on  the  signal  variance.  Specifically, 
consider  the  first-order  network  drawn  in  Fig.  5  with  the  two 
noise  sources  ei(«)  and  etin).  From  the  model  for  multiplier 


•t  Is) 


Fig.  5.  Noisy  first-order  filter  (floating  point). 


roundoff  noise,  the  noise  source  Ci(n)  is  given  by 


din)  =  ayin  -  1)«. 

(33) 

and  the  noise  source  etin)  is  given  by 

*i  i")  =  gin)£„. 

(34) 

The  analysis  above  in  which  we  neglected  second-order  terms 
corresponds  in  this  case  to  evaluating  the  variance  of  ei(n)  and 
ej(»)  by  using  the  mean-square  values  for  yin—  1)  and  g in) 
that  would  result  if  no  roundoff  noise  were  present.  Therefore, 
if  we  assume  that  x(n)  is  a  zero-mean  white-noise  input,  with 
variance  oS.  then  the  variances  of  «i(n)  and  etin)  are,  respec¬ 
tively, 

a,,*  =  aV#2v*(w  —  1)  =  - -  (35) 

1  —  ar 

=  or{y(w)  =  a*V*2  - - -  (36) 

1  —  a2 

where  the  bar  denotes  expected  value.  Then,  since  ei(«)  and 
etin)  are  independent,  because  t„  and  £„  are  inde|>endent,  the 
output  noise  variance  is 

1  +  a'-  1  +  a' 

”o  -  — - —  =  *.*»,*  - - -  (37a) 

(1  —  a*)2  1  —  a2 

where  we  have  assumed  again  that  <r,2  and  <T{2  arc  equal.  The 
output  noise-to-signal  ratio  is 


a  o' 

17„2 


1  +  a2 
1  —  a- 


(37b) 


We  ran  analyze  the  effect  of  roundoff  noise  in  the  second-order 
filler  in  a  similar  manner.  In  Fig.  6  is  shown  the  network  for  a 
second-order  filter  with  roundoff  noise  sources  included.  Note 
that  since  noise  sources  must  lie  included  due  to  addition,  two 
summers  are  included  to  add  the  three  variables  in  the  feed¬ 
back  loop.  The  noise  sources  ej(n)  and  cpn)  represent  the  noise 
due  to  the  multiplies  and  the  noise  sources  edn)  and  *•(«)  rep¬ 
resent  the  noise  due  to  the  additions.  With  assumptions  simi¬ 
lar  to  those  above  in  w  hich  we  neglected  second-order  terms, 
we  write  that 


tiin)  =  yin)n(n) 

<’••!(»)  ”  [.'•(«)  -  *(»)]«»(«) 
e>(n)  =  2r  cos  8y(n  —  1  )«3(h) 
e,in)  =  -  r2y(w  -  2)h(h) 


S4 


i 


k 


(38) 


■"•Mann 
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where  «>. «j,  nml  «,  are  indc|H*iidcnt  random  variables  with 
e<|nal  variance  <r,1.  If  as  l>cfore,  x(n)  is  assumed  to  be  a  while 
random  process  with  variance  a,',  then  the  output  noise-to- 
ti|{nal  ratio  for  the  second-order  case  is 

o.7  f  /  r*  cos’  0\T 

— “  “  •+■  6'^3r4  +  12r*  tos’0  —  16  — — ~t)J  (*^) 

where 

c  1  +  ft  /  I  \ 

1  -  r1  \r4  +  1  -  4r5 cos1 6  +  2r’/'  (40) 

bor  the  high  gain  case,  it  is  (Kissible  to  compare  fixed-point 
and  ll< >ii ti n vt - p->i n t  arithmetic  by  approximating  the  expres¬ 
sions  for  the  imisc-to-simial  ratio.  For  the  first-order  case, 
with  it-l -4.  and  4«1.  the  remit  U7I.)  for  the  first-order 
filter,  can  Ik-  approximated  as 


o.7  1  I 


(41) 


Similarly  for  the  second  order  filter 


/.l  +  4  cos’  0\ 
\  44  sin*  6  ) 


(42) 


where  in  (41)  and  (42)  we  have  taken  <r.*  —  (1/3)2“^. 

l  or  lixed-|Kiint  arithmetic  we  recall  that  fora  white-noise 
input  the  noise-  to-signal  ratio  India  veil  as  1/4’ and  fora  sinus¬ 
oidal  input  as  I  4.  Comparison  of  (41)  and  (42)  with  (14)  and 
'  lb)  and  (IK)  and  (20)  indicates  a  slightly  larger  uoisc-to-sig- 
mil  ratio  for  llontiug-|>oitil  arithmetic  as  compared  with 
fixcd-|K>iut  arithmetic  with  a  sinusoidal  input  of  known  fre- 
tpiency  Imt  a  significantly  smaller  noise- to-signal  ratio  for 
lloutiiiK-iMiiut  arithmetic  as  compared  with  fixed-|Miint  arith¬ 
metic  with  a  white-noise  input  It  is  inijiorlant  to  keep  in 
mind,  however,  that  the  uoiso-io-signal  ratios  for  the  fixed- 
|M.iin  filters  were  computed  on  the  basis  that  the  input  signal 
was  as  large  as  posmhle  If  the  inpui  signal  level  decreases,  the 
noise- to-sigual  ratio  will  increase  since  the  output  noise  vari- 
-nice  is  independent  of  the  input  signal  level.  For  n<.ating-|»oiiii 
aritli uietic,  on  the  other  hand,  the  output  noise  variance  is 
pro|Hirtional  to  the  output  signal  variance  and  as  the  input 
level  is  scaled  up  or  down  so  is  the  roundoff  noise.  Il  is  also 
ini|sirianl  to  note  that  the  comparison  just  discussed  assumes 
that  the  11oHiing-|Niiut  mantissa  is  equal  in  length  to  the 
intire  fixed  |>oini  word,  and  does  not  account  for  the  extra 
bits  needed  for  the  characteristic  1  he  authors  [6]  have  previ¬ 


ously  compared  fixed-  and  lloating-point  filters  on  the  basis  of 
equal  total  word  length.  However,  in  completing  such  a 
comparison  one  must  take  account  of  the  large  difference  in 
hardware  complexity  between  implementing  lloaling-|M)int 
arithmetic,  and  adding  a  few  bits  to  a  fixed -point  arithmetic 
element. 

OpiH'idicim  |0)  has  proposed  a  realization  of  recursive 
digital  filters  using  block  floating-point  arithmetic.  Here  the 
input  and  filter  states  (i.e„  the  inputs  to  the  delay  registers) 
are  jointly  normalized  before  the  multiplications  and  addi 
lions  are  performed  in  fixed-point  arithmetic.  The  scale  factor 
(or  exponent)  obtained  during  the  normalization  is  then  ap¬ 
plied  to  the  final  output  to  obtain  a  fixed -point  output.  The 
roundoff  noise*  properties  of  such  a  realization  were  studied, 
and  the  noise-to-signal  ratio  was  found  to  lie  between  that 
for  fixed  and  floating  |M)int. 

h.  Zero-Input  Limit  Cycle  Behavior  of  Digital  hitters  for 
Fixed- Point  .1  ritlimelic 

In  the  preceding  discussion  the  effect  of  arithmcti  round- 
olf  was  modelled  as  an  additive  white-noise  source,  tincorro- 
Inicd  with  the  data.  Justification  of  this  model  assumes  that 
from  iteration  to  iteration,  the  input  can  be  expected  to  pass 
through  several  quantization  levels.  Consequently,  this  model 
is  applied  primarily  when  the  input  signal  has  a  complicated 
behavior  and  cannot  be  expected  to  be  valid  in  general.  For 
example,  consider  a  lirst-order  filter  for  which  the  difference 
equation  is 

V.  =  oy„-i  +  x„  (43) 

and  for  which  the  register  length  for  the  data  is  4  bits  and  the 
coefficient  a  is  0.5.  If  the  input  x„  is  7/8  and  if  rounding  is 
applied  after  the  arithmetic  then  on  successive  iterations  of 
the  filter,  the  output  will  he: 

)'o  =  7/8 

)'t  =  1/2 

V:  ■»  1/4 

>'»  =  1/8 

y*  =  1/8,  for  n  >  4. 

Thus  due  to  rounding,  the  output  reaches  a  steady-state  non¬ 
zero  value  ami  since  the  ideal  steady-state  output  is  zero,  this 
nonzero  value  represents  roundoff  error.  Clearly  this  kind  of 
roundoff  error  cannot  Ik*  modelled  as  white  noise,  but  in  fact 
represents  a  limit  cycle  due  to  the  nonlinearity  corresponding 
to  the  quantizer  which  implements  the  rounding.  Limit  cycle 
behavior  of  this  type  was  first  noted  by  Blackman  ( 10]  who 
referred  to  the  amplitude  intervals  within  which  these  limit 
cycles  are  conlincd  as  “deadbands.”  Blackman  considered 
only  lirst-order  limit  cycles  corresponding  to  a  ilc  behavior  in 
tile  dead  hand.  More  generally,  Jackson  |ll)  has  considered 
limit  cycle  behavior  in  first-  and  second-order  filter  sections 
with  an  analysis  based  on  the  location  of  the  “effective”  poles 
in  ilu*  filter  due  to  roundoff.  Following  the  approach  pre¬ 
sented  hv  Jackson,  consider  a  first-order  filter  with  a  differ 
cnee  equation  of  the  form  of  (43).  Due  to  the  register  length 
constraint,  the  product  av„_i  must  he  rounded.  Let  (•)' 
denote  the  iqKration  of  rounding.  If  the  register  length  is 
(4  +  1)  hits  and  if  data  are  represented  as  fractions  then 

|  (ay„-,)'  -  ay„_,  <  (J)2-‘.  (44) 

II  }*.-i  is  such  that  (ay,,_i)'|  =  |v„_i  then  the  magnitude  of 
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the  effective  value  of  the  coefficient  is  unity  correspond! ng  to 
the  pole  of  the  filter  being  on  the  unit  circle.  The  range  of 
values  for  which  this  condition  is  met  is 


y-A  ~  kv„-,|  <(})2-k 

(45) 

(0.5) 2-" 

>■-  £i-i„r 

(46) 

This  range  of  values  is  referred  to  as  the  deadband.  Due  to 
rounding,  of  course,  values  within  the  dead  hand  must  he  in 
steps  of  2  'b.  For  the  first-order  filter,  when  the  filter  state 
falls  within  this  range  and  the  input  is  zero,  the  effective  pole 
is  on  the  unit  circle  and  the  filter  will  support  a  limit  cycle 
behavior.  If  the  coefficient  a  is  positive,  as  in  the  above  ex¬ 
ample,  the  limit  cycle  response  is  dc,  i.e.,  has  constant  mag¬ 
nitude  and  sign.  For  a  negative  the  limit  cycle  behavior  has 
constant  magnitude  hut  alternating  sign. 

For  a  seemd-order  filter  there  is  a  larger  variety  of  modes 
of  limit  cycle  behavior.  In  particular,  consider  the  second- 
order  difference  equation 

y„  =  xn  ~  /3|V„_i  -  /32yn_5.  (47) 


With  di!<+4/32  the  filter  poles  occur  as  a  complex  conjugate 
pair  and  with  /3j=l  the  poles  occur  on  the  unit  circle.  The 
approach  proposed  by  Jackson  for  examining  the  limit  cycle 
behavior  of  the  second-order  filter  corresponds  to  considering 
the  filter  behavior  when  the  effect  of  rounding  places  the  effec¬ 
tive  poles  of  the  filter  on  the  unit  circle.  With  zero  input  the 
effective  poles  will  he  on  the  unit  circle  if 

I  v„— 2 1  -  I  (&v„_..)  |  <  •}  2-‘  (48) 


or 


Vn-2  I  < 


(0.5) 2‘ A 

M/M  ‘ 


(49) 


Thus  if  the  output  falls  within  this  range  the  effective  value  of 
0i  is  unity  so  that  the  effective  poles  are  on  the  unit  circle. 
With  the  effective  value  of  0»  as  unity,  the  effective  value  of 
0i  controls  the  oscillation  frequency. 

A  second  mode  of  limit  cycle  behavior  occurs  in  second- 
order  filters  when  the  effect  of  rounding  is  to  place  an  effective 
pole  at  s=±1  As  shown  by  Jackson,  the  deadband  corre¬ 
sponding  to  this  mode  is  for  values  less  than  oi  equal  to 
1  (1  — |/Ji|  +0i)  in  steps  of  integer  multiples  of  2  (l. 

While  this  approach  is  somewhat  heuristic,  Jackson  has 
found  that  these  bounds  are  consistent  with  experimental 
results  and  hence  he  has  hypothesized  that  they  represent 
necessary  and  sufficient  conditions.  These  bounds  for  second- 
order  filters  are  summarized  in  Fig.  7,  showing  different  dead¬ 
band  subregions  in  the  0\,  0i  plane.  The  number  within  an 
area  in  the  0i,  0:  plane  represents  the  maximum  magnitude  of 
the  limit  cycle  in  multiples  of  2~b  and  the  cross  hatched  region 
represents  the  region  for  which  no  limit  cycles  can  occur. 

Recently,  Parker  and  Hess  1 1 2 ]  have  studied  the  limit 
cycle  problem  further,  and  found  that  these  bounds  are  ap¬ 
proximately  correct  and  sufficient,  hut  not  necessary.  In 
other  words,  there  exist  some  limit  cycles  outside  the  regions 
specified  by  Fig.  7. 

In  addition  to  the  above  classes  of  limit  cycles,  a  more 
severe  tvpe  of  limit  cycle  can  occur  due  to  overflow  in  filters 
implemented  using  one’s-coniplenient  or  t wo’s-complemcnt 


arithmetic.  These  limit  cycles  have  been  referred  to  as  over- 
llovv  oscillations  [  1 3 ]  and  can  be  avoided  by  using  saturation 
arithmetic. 

F.  Effects  of  Parameter  Quantization  in  Digital  Filters 

In  the  preceding  sections  we  focussed  on  the  effects  of 
arithmetic  roundoff  in  digital  filters.  Another  consequence  of 
the  requirement  of  finite  register  length  is  that  the  filter  coeffi¬ 
cients  cannot  he  specified  exactly.  Classical  design  procedures 
generally  lead  to  filter  coefficients  with  arbitrary  accuracy  and 
the  implementation  of  the  filter  then  requires  that  the  coeffi¬ 
cients  he  modified  to  lit  the  available  register  length.  For  hard¬ 
ware  realizations  of  digital  filters  it  is,  of  course,  desirable  to 
keep  the  register  length  as  small  as  possible. 

One  common  approach  to  the  problem  of  parameter  quan¬ 
tization  is  the  use  of  filter  configurations  or  structures  which 
in  some  sense  are  least  sensitive  to  inaccuracies  in  the  param¬ 
eters.  One  of  the  difficulties  in  evaluating  the  sensitivity  of 
filter  structures  is  the  choice  of  a  meaningful  measure  of  the 
sensitivity.  Most  commonly,  the  sensitivity  of  the  filter  is  tied 
to  the  movement  of  the  poles  of  the  filter.  For  this  choice 
Kaiser  [  1 4 ]  has  shown  that  for  a  filter  with  clustered  poles  a 
cascade  or  parallel  combination  of  first-  and  second-order  sec¬ 
tions  provides  more  accuracy  in  the  pole  positions  than  a 
direct  form  realization.  This  is  basically  a  consequence  of  the 
fact  that  for  a  polynomial  whose  roots  are  clustered,  the  sensi¬ 
tivity  of  the  roots  to  changes  in  the  polynomial  coefficients 
increases  as  the  order  of  the  polynomial  increases.  Thus  the 
roots  can  he  more  accurately  controlled  if  the  polynomial  is 
factored  into  first-  and  second-order  factors. 

F.ven  within  the  choice  of  first-  and  second-order  sections 
some  flexibility  remains.  For  a  direct  form  implementation  of 
a  pole  pair  as  shown  in  Fig.  8(a)  the  coefficients  are  —  r1  and 
2 r  cos 0.  For  a  given  quantization  on  the  coefficients  the  poles 
must  lie  on  a  grid  in  the  z  plane  defined  by  the  intersection  of 
concentric  circles,  cnrrespondi»g  to  quantization  of  r*  and 
vertical  lines,  corresponding  piantization  of  2r  cos  0  Such 
a  grid  is  illustrated  in  Fig.  8(lq.  An  alternative  realization  of 
a  pole  pair  is  the  coupled  form  proposed  by  Rader  and  ( ioltl 
[l5],  as  shown  in  Fig.  9(a).  In  this  case  the  coefficients  are 
r  cos  0  and  r  sin  0  and  consequently  the  poles  must  lie  on  a 
rectangular  grid  as  illustrated  in  Fig.  9(h).  We  note,  for  ex¬ 
ample,  that  fora  given  coefficient  word  length  the  direct  form 
permits  more  accurate  placement  of  poles  with  r  close  to  unity 
and  0  large  while  the  coupled  form  is  more  advantageous  for  0 
small.  There  arc,  in  theory,  many  other  structures  in  addi¬ 
tion  to  the  direct  and  coupled  forms  for  implementing  pole 
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Mg.  H.  (a)  Direct  form  Implementation  of  a  pole  pair. 
(1>)  Grid  of  allowanle  pole  positions—  direcl  form. 


pairs  although  they  are  the  most  commonly  considered  [  16]. 
Different  structures,  of  course,  imply  different  grids  in  the  s 
plane  and  generally  it  is  advantageous  to  choose  a  structure 
for  which  the  grid  is  dense  in  the  region  of  the  s  plane  where 
the  poles  are  to  lie  located. 

With  a  given  choice  of  structure  there  remains  the  ques¬ 
tion  as  to  how  the  pole  locations  on  the  grid  should  lie  chosen. 
A  common  procedure  is  to  truncate  or  round  the  ideal  coeffi¬ 
cients.  An  alternative  used  by  Avenhaus  and  Schussler  [ 1 7 ] 
and  also  by  Steiglitz  [  1 8 ]  is  to  search  over  the  grid  in  the  vi¬ 
cinity  of  the  ideal  pole  locations  to  select  a  grid  point  which 
locally  minimizes  the  maximum  error  in  the  filter  frequency 
response.  As  an  alternative  to  the  use  of  cascade  or  parallel 
connections  of  first-  and  second-order  sections,  more  general 
filter  structures  can  be  considered.  Digital  wave  filters,  as 
proposed  by  Fettxveis  [19]  and  investigated  by  Bingham  [20 ] 
and  by  Crochiero  [21  ],  appear  to  have  much  less  sensitivity 
to  parameter  inaccuracies  than  the  cascade  form. 

It  would,  of  course,  be  desirable  to  incorporate  the  con¬ 
straint  of  quantized  coefficients  into  the  design  of  digital  fil¬ 
ters.  For  nonrecursive  filters,  algorithmic  design  falls  within 
the  framework  of  integer  linear  programming.  For  recursive 
filters,  however,  the  equations  become  nonlinear.  In  general, 
the  development  of  design  procedures  with  quantized  coeffi¬ 
cients  remains  an  important  area  of  research. 

IV.  Effects  of  Arithmetic  Roundoff  in  the  FFT 
.■I.  Introduction 

The  FFT  algorithm  [22]  for  computing  the  discrete  Fou¬ 
rier  transform  ( DFT)  plays  a  central  role  in  many  signal  pro¬ 
cessing  applications  [23],  As  with  the  implementation  of 
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(a) 


(b) 


Fig.  9.  (si)  Coupled  form  iuiplemenlalion  of  a  pole  pair. 

(b)  ('.rid  of  allowable  pole  posit  ions— coupled  form. 

digital  filters,  it  is  important  to, understand  the  effect  of  finite 
register  length  arithmetic  on  the  performance  of  the  algo¬ 
rithm. 

There  are  many  forms  of  the  FFT  algorithm  and  the  de¬ 
tailed  effects  of  quantization  will  differ  depending  on  the 
form  used.  The  most  commonly  used  forms  of  the  algorithm 
are  the  radix-2  forms  for  whicl  the  size  of  the  transform  com¬ 
puted  is  an  integer  power  of  tv/o.  For  the  most  part,  the  dis¬ 
cussion  below  is  phrased  in  ter  lis  of  a  particular  form  of  the 
radix-2  FFT,  commonly  referred  to  as  the  decimation  in  time 
form  of  the  algorithm;  the  results  however  are  applicable  with 
only  minor  modification  to  the  decimation  in  frequency  form. 
We  feel  that  most  of  the  ideas  employed  in  the  error  analysis 
of  the  radix-2  forms  of  the  algorithm  can  be  utilized  in  other 
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forms  such  as  mixed  radix,  etc. 

Our  approach  in  analyzing  noise  in  the  FFT  is  basically 
statistical.  In  most  cases,  the  predictions  of  the  models  are 
supported  with  experimental  data  (from  Weinstein  [24],  un¬ 
less  otherwise  stated).  For  floating  and  block  floating  point 
arithmetic,  in  order  to  simplify  the  analysis  and  obtain  con¬ 
crete  results,  it  is  convenient  to  assume  a  simple,  white-noise 
model  for  the  signal  being  transformed.  Discussion  of  how  the 
results  might  be  expected  to  change  for  other  types  of  signals 
is  included,  as  are  experimental  noise  measurements  on  FF I 's 
of  nonwhite  signals. 

B.  The  FFT  Algorithm 

The  FFT  algorithm  is  directed  toward  computing  the 
DFT  of  a  finite  duration  sequence /(«),  defined  as 


/ 
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l*'ig.  10.  KKT  How  grapli.  A'  =  H. 


,V-I 

p(k)  ~  X  ./(»)  if  =  c~i{triX).  (50) 

n-(l 

A  flow  chart  depicting  the  FFT  algorithm  for  N=S  =  23  is 
shown  in  Fig.  10,  A  specific  rlcciination  ill  time  algorithm  is 
depicted.  (An  implementation  of  this  particular  form  of  the 
algorithm  was  used  for  the  reported  experimental  work.) 
Some  key  aspects  of  this  diagram,  which  are  common  to  all 
standard  radi.\-2  algorithms,  are  as  follows.  The  DFT  is  com¬ 
puted  ill  F‘=logs  A  stages.  At  each  stage,  the  algorithm  passes 
through  the  entire  array  of  A'  complex  numbers,  two  at  a 
time,  generating  a  new  N  number  array.  'Flic  eth  array  con¬ 
tains  the  desired  DFT.  The  basic  numerical  computation 
operates  on  a  pair  of  numbers  in  the  (»/  +  l)tli  array.  This 
computation,  referred  to  as  a  “butterfly”  is 

A’„.h(0  =  A '„,(/)  H-  IF.V„,(i) 

A’mnO)  =  A’m(»)  -  \VXm(j).  (51) 

Here,  A'„,(t)  and  A’„,(/)  represent  a  pair  of  numbers  in  the  wth 
array,  and  IT’  is  some  appropriate  integer  power  of  IF,  that  is 

IF  =  JF»-  =  c->Sr"/v.  (52) 

The  form  of  the  butterfly  computation  is  actually  somewhat 
different  for  a  decimation  in  frequency  algorithm,  where  the 
computation  is 

A',„h(i)  =  A'm(t)  +  Xm(j) 

A’.,(  ,(i)  =  l-Ym(i)  -  A'm(i)]lF.  (53) 

At  each  stage,  N/2  separate  butterfly  computations  are  car¬ 
ried  out  to  produce  the  next  array.  The  integer  p  varies  with 
i,  j,  and  m  in  a  manner  which  depends  on  the  specific  form  of 
the  FFT  algorithm  that  is  used.  Fortunately,  our  analysis  is 
not  tied  to  the  specific  way  in  which  p  varies.  Also,  the  specific 
relationship  between  i,  j,  and  m,  which  determines  how  we 
index  through  the  with  array,  is  not  important  for  the  anal¬ 
ysis.  The  details  of  the  analysis  for  decimation  in  time  and 
decimation  in  frequency  differ  somewhat  due  to  the  different 
butterfly  forms,  but  the  basic  results  for  the  dependence  of 
noise-to-signal  ratio  on  N  do  not  change  significantly.  In  our 
analysis  we  will  assume  a  butterfly  of  the  form  (51),  corre¬ 
sponding  to  decimation  in  time. 


C.  FFT  Roundoff  Noise  with  Fixed-Point  Arithmetic 

We  will  model  the  roundoff  noise  by  associating  an  inde¬ 
pendent  wdiite-noise  generator  with  each  multiplier.  This 
means  that  a  noise  source  feeds  into  each  node  of  the  signal 
(low  graph  of  Fig.  10  (excluding  ihe  initial  array  of  nodes, 
since  we  are  not  considering  A/D  noise  here).  Since  we  are 
dealing  with  complex  multiplications,  these  elemental  noise 
sources  are  complex.  Defining  the  complex  variance  o/r  as  the 
expected  squared  magnitude  of  such  a  noise  source,  we  have 

2-26 

air  =  4  — —  (54) 

12 

where  it  is  assumed  that  each  of  the  four  real  multiplications 
used  to  perforin  the  complex  multiplication  is  rounded  sepa¬ 
rately.  In  F'ig.  10,  3X8  =  24  such  noise  sources  must  be  in¬ 
serted.  To  add  the  effects  of  each  of  the  noise  sources  in  eval¬ 
uating  the  total  roundoff  noise  in  the  output,  we  note  that  the 
transmission  function  from  any  node  in  the  flow  graph  to  any 
other  connected  node  is  multiplication  by  a  complex  constant 
of  unity  magnitude.  Since  we  assume  that  all  noise  sources  are 
uncorrclatcd,  the  noise  variance  at  any  output  node  is  equal  to 
an1  times  the  number  of  noise  sources  that  propagate  to  that 
node.  The  general  result  which  is  easily  verified  for  the  case 
N=  8  by  inspection  of  F'ig.  10  is  that  (A-l)  noise  sources 
propagate  to  each  output  node  so  that  the  output  noise  vari¬ 
ance  <rp*  is  given  by 

air  =  (Ar  -  1  )a,r- 

wliich  for  large  N  we  take  as 

aff1  ~  Na„2.  (55) 

According  to  this  result,  the  variance  of  the  output  noise  is 
proportional  to  N,  the  number  of  points  transformed.  The 
effect  of  doubling  N,  or  adding  another  stage  in  the  FFT,  is  to 
double  the  output  noise  variance.  Using  the  assumptions  we 
have  made  thus  far  about  the  noise  generators  in  the  FFT  (all 
uncorrclatcd  with  equal  variances),  the  output  noise  is  white, 
i.c.,  the  N  noise  samples  E{k)  are  mutually  uncorrelateci,  with 
iiidciiciidcnt  real  and  imaginary  parts.  This  follows  from  the 
fact  that  the  output  of  any  butterfly  is  white  (two  outputs  uii- 
corrt!  ited  with  equal  variance,  real  and  imaginary  parts 
uncorrclatcd)  if  the  input  is  white.  Since  the  noise  sources  in 
our  system  are  white,  and  all  connected  to  the  output  via 
some  combination  of  butterfly  computations,  the  output  noise 
must  also  be  white. 

In  order  to  simplify  the  analysis  leading  to  (55),  we  have 
neglected  some  details.  First,  we  have  associated  equal  vari¬ 
ance  noise  sources  with  all  multipliers,  including  where  IF=1 
and  j.  In  many  programmed  FFT’s  these  multiplications  are 
performed  noiselessly.  If  we  assume  in  the  analysis  that  these 
multiplications  are  noiseless,  the  output  noise  variance  will  no 
longer  be  uniform  over  the  output  array.  For  example,  the 
zeroth  output  point  would  be  noiseless.  The  average  variance 
over  the  output  array  will  be  somewhat  lower  than  the  result 
in  (55),  but  will  retain  a  linear  dependence  on  N.  Second,  the 
assumption  that  all  noise  sources  are  uncorrelated  is  contra¬ 
dicted  by  the  fact  that  the  two  noise  sources  associated  with  a 
given  butterfly  arc  negatives  of  each  other,  and  therefore  com¬ 
pletely  correlated.  This  docs  not  affect  the  result  for  output 
noise  variance,  since  the  two  outputs  of  a  butterfly  connect  to 
a  disjoint  set  of  output  points.  However,  it  implies  that  the 
output  noise  samples  E{k)  are  somewhat  correlated.  These 
details  are  worth  mentioning,  but  not  worth  analyzing  here  at 
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length,  because  they  cloud  the  essential  ideas  of  the  analysis, 
are  quite  program-dependent,  and  do  not  change  the  essential 
character  of  the  dependence  of  mean-squared  output  noise 
on  N. 

In  implementing  the  FFT  with  fixed-point  arithmetic  we 
must  insure  against  overflow.  From  (51)  it  follows  that 

max  [  |  A'm(i)  |  ,  |  X ,„(;)  |  ] 

<  max  [  |  A'mt-i(i)  |  j  |  |  ]  (56) 

and  also  that 

max  [|  A'm ,  ,(;)  | ,  |  A'm,  i(j)  |  ] 

<  2  max  [  |  A'm(t)  |  ,  |  Xm(j)  |  ].  (57) 

Equation  (56)  implies  that  the  maximum  modulus  is  non¬ 
decreasing  from  stage  to  stage  so  that,  if  the  magnitude  of  the 
output  of  the  FFT  is  less  than  unity  then  the  magnitude  of 
the  points  in  each  array  must  be  less  than  unity,1  i.e.,  there 
will  be  no  overflow  in  any  of  the  arrays. 

In  order  to  express  this  constraint  as  a  bound  on  the  input 
sequence,  we  note  that  the  maximum  possible  output  can  be 
expressed  in  terms  of  the  maximum  input  as 

.v-i 

Af(  k)  | . *  <  |  ,r(»)  I,,,,,  £  j  IF"*  |  =  .V  |  x(n)  [,„»„.  (58) 

ri-0 

Thus  hounding  the  input  sequence  so  that 

|  *(>«)  |  <  1/A'  (59) 

will  prevent  overflow.  To  obtain  an  explicit  expression 
for  output  signal  variance,  we  assume  .v(w)  white,  with 
real  and  imaginary  parts  each  uniformly  distributed  in 
(  —  1/ \/2 A',  1/%/2-V).  Then  we  have 

ir.\-  =  |  X(W  =  AV,*  =  .v  1707)1 2  =  -T  •  (60) 

Combining  this  with  (55)  yields 

=  3.V-V/I2.  (61) 

Os' 

The  assumption  of  white  input  signal  is  not  critical  here.  For 
example,  if  a  complex  sinusoid  x(n)  =  (1/ N)  cx\ij(2wko>i/N-\-<t>) 
had  been  selected  Oh}/os2  would  still  be  proportional  to  A’2, 
which  is  the  essential  point  of  (61). 

Equation  (57)  suggests  an  alternative  procedure  for  pre¬ 
venting  overflow.  Since  the  maximum  modulus  increases  by 
no  more  than  a  factor  of  two  from  stage  to  stage  we  can  pre¬ 
vent  overflow  by  requiring  that  |.v(«)|  <1  and  incorporating 
an  attenuation  factor  of  1/2  at  each  stage.  Using  this  step-by- 
step  scaling,  the  attainable  output  signal  level  (for  white 
input)  is  the  same  as  in  (60)  since  the  output  signal  level  does 
not  depend  on  where  the  scaling  is  done,  but  only  on  how 
much  overall  scaling  is  done.  However,  the  output  noise  level 
will  bi“  much  less  than  in  (55)  since  the  noise  introduced  at 
early  stages  of  the  FFT  will  be  attenuated  by  the  scaling 
which  takes  place  at  the  later  array.  Quantitatively  for  N  =  2" 

1  Act  unity  mi  should  discuss  overflow  ill  terms  of  the  real  and  imiiRl- 
nary  parts  ol  the  data,  rather  than  the  magnitude.  However,  x  <1 
implies  that  Ui  (x)  <1  and  !m  (x)  |  <1,  and  only  a  slight  Increase  in 
allowable  signal  level  is  achieved  by  scaling  on  the  basis  of  !<c  and  Im 
parts. 
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of2  =  <r/i'2  S  7  (62) 

*-l  k 

where  <r/e2  represents  the  roundoff  noise  introduced  due  to 
multiplication  by  IF  and  scaling  and  will  consequently  be 
slightly  higher  than  o/}.  In  particular,  if  we  assume  that  the 
scaling  is  accomplished  with  rounding,  it  can  he  shown 

5 

o if*  =  -  l~2b.  (63) 

6 

For  large  N,  (62)  is  approximately 

of}  =  2<t;j.2  (64) 

and  thus  is  much  less  than  the  noise  variance  resulting  when 
all  of  the  scaling  is  carried  out  on  the  input  data. 

Now,  wc  can  combine  (64)  with  (60)  to  obtain  the  output 
noise-to-signal  ratio  for  the  case  of  step-by-step  scaling  and 
white  input.  We  obtain 

Oh} 

- =  6.Y<r/i-2  =  (5  A”)  2— 26  (65) 

ox2 

a  result  proportional  to  N,  rather  than  to  N1.  An  interpreta¬ 
tion  of  (65)  is  that  the  rms  output  noise-to-signal  ratio  in¬ 
creases  as  N,  or  by  half  a  bit  per  stage.  This  result  was  first 
obtained  by  Welch  [25].  It  is  important  to  note  that  the  as¬ 
sumption  of  white  signal  is  not  essential  in  the  analysis.  The 
basic  result  of  half-a-bit-per-stage  increase  holds  for  a  broad 
class  of  signals,  with  only  the  constant  multiplier  in  (65) 
being  signal-dependent.  In  particular,  for  a  general  input  with 
settling  at  each  array,  the  output  variance  is  related  to  the 
variance  of  the  input  array  by 

1  1 _ 

a. y2  —  —  o}  =  —  |  x(ii) ;  -  (66) 

so  that 

-  X2~2b 

oy}  3 

—  = -  (67) 

OX2  o2 

where,  to  reduce  noise-to-signal  ratio,  we  would  like  to  make 
o }  as  large  as  possible  but  are  limited  by  the  constraint 
|.v(»»)|  <1.  The  result  (67)  has  been  verified  experimentally 
for  both  wide-hand  and  narrow-hand  signals  [24],  [25]. 

We  should  also  note  that  the  dominant  factor  causing  the 
increase  of  opr  ox2  with  N  is  the  decrease  in  signal  level 
(required  by  the  overflow  constraint)  as  we  pass  front  stage  to 
stage.  According  to  (63)  and  (64),  very  little  noise  (only  a  bit 
or  two)  is  present  in  the  final  array.  Most  of  the  noise  has 
been  shifted  off  by  the  sealings.  However,  the  mean-squared 
signal  level  lias  decreased  by  <t  factor  of  1  / N  from  its  initial 
value,  due  to  the  scalings.  Our  output  consists  not  of  the  DFT 
defined  by  (50)  but  of  l/A7  times  this  DFT. 

We  have  assumed  straight  fixed  point  computation  in  this 
section,  i.c.,  only  preset  attenuations  were  allowed,  and  we 
were  not  permitted  to  rescale  on  the  basis  of  an  overflow  test. 
Clearly,  if  the  hardware  or  programming  facility  arc  such  that 
straight  fixed  point  must  lie  used,  wc  should,  if  possible,  incor¬ 
porate  attenuators  of  1/2  at  each  array  rather  than  using  a 
largt  attenuation  of  the  input  array. 
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log.  1 1.  Experimental  und  theoretical  noise-to-signal 
ratios  for  block  floating-point  KHT. 


A  third  approach  to  avoiding  overllow  is  the  use  of  block 
lloating  point.  In  this  procedure  the  original  array  is  normal¬ 
ized  to  the  far  left  of  the  computer  word,  with  the  restriction 
that  |  x(h  <1;  the  computation  proceeds  in  a  fixed  point 
manner,  except  that  after  every  addition  there  is  an  overllow 
test;  if  overllow  is  detected,  the  entire  array  is  shifted  right  1 
hit  and  the  computation  continues.  The  number  of  necessary 
shifts  are  counted  to  determine  a  scale  factor  or  exponent  for 
the  entire  final  array.  The  output  noisc-to-signal  ratio  depends 
strongly  on  how  many  overflows  occur,  and  at  what  stages  of 
the  FFT  they  occur.  The  positions  and  timing  of  overflows 
are  determined  by  the  signal  being  transformed,  and  thus,  in 
order  to  analyze  noise-to-signal  ratio  in  block  floating  FFT, 
one  needs  to  know  the  signal  statistics.  This  is  in  contrast  to 
the  fixed  point  analysis  above,  where  it  was  not  necessary  to 
assume  specific  signal  statistics. 

The  necessary  number  of  right  shifts  of  the  array  is  related 
to  the  peakiness  of  the  DFT  of  the  signal  being  transformed. 
If  the  constant  signal  ,v(»)  =  1  or  the  single-frequency  input 
.v(n)=  ux\>  j(2ir  /  N)k0n  is  transformed,  the  output  (with  k„  an 
integer)  will  consist  of  a  single  nonzero  point  and  (for  N  =  2’)  v 
scalings  of  the  array  will  be  necessary,  one  for  each  stage. 

A  reasonable  case  to  examine  is  the  case  of  a  white  input 
signal;  the  DFT  of  a  white  signal  is  white,  and  one  might 
expect  (since  the  spectral  energy  is  spread)  that  scalings  at  all 
stages  would  not  he  necessary,  and  a  noisc-to-signal  ratio 
advantage  over  fixed  point  would  be  gained.  This  problem 
can  be  analyzed  t!  corctically  [24]  but  the  analysis  is  quite 
involved  and  will  be  omitted.  Instead,  we  will  present  some 
experimental  results. 

In  Fig.  11  experimentally  measured  values  of  output 
noise-to-signal  ratio  arc  presented  for  block  floating  FF'T’s  of 
white  inputs,  using  rounded  arithmetic.  The  quantity  plotted 
is  [off1  2"2VI2)1'2,  the  r ms  noise-to-signal  ratio.  For  com¬ 
parison.  a  theoretical  curve  representing  fixed  point  noisc-to- 
signal  ratio  (for  rounded  arithmetic)  is  also  shown.  We  see 
that  for  white  input  block  lloating  point  provides  some  ad¬ 
vantages  over  fixed  point,  especially  for  the  larger  transforms. 
For  Ar=2048,  the  rms  noise-to-signal  ratio  for  block  lloating 
point  is  about  1  8  that  of  fixed  point,  representing  a  3-hit 
improvement. 

An  experimental  investigation  was  used  to  examine  how 


Fig.  12.  Noisy  butterfly  computation  (floating  point). 


the  results  for  block  lloating  point  change,  when  truncation 
rather  than  rounding  is  used.  The  results  of  this  experiment 
are  also  shown  in  Fig.  11.  Noise-to-signal  ratios  are  generally 
a  bit  or  two  worse  than  for  rounding.  The  rate  of  increase  of 
noise-to-signal  ratio  with  N  seems  to  be  about  the  same  as  for 
rounding. 

D.  FFT  Roundoff  Noise  with  Floating-Point  Arithmetic 

The  effect  of  arithmetic  roundoff  with  floating-point  arith¬ 
metic  has  been  analyzed  theoretically  and  experimentally  by 
Gentleman  and  Sande  [26],  by  Weinstein  [27],  and  by 
Kancko  and  Liu  [28].  As  with  the  statistical  analysis  of 
roundoff  errors  with  fixed-point  arithmetic,  noise  is  intro¬ 
duced  due  to  each  butterfly  computation.  As  with  floating¬ 
point  errors  in  digital  filters,  we  neglect  second-order  error 
terms  so  that  noise  sources  are  introduced  after  each  multi¬ 
plication  and  addition  that  are  assumed  to  be  white  but  for 
which  the  variance  is  proportional  to  the  variance  of  the  sig¬ 
nal  at  that  node.  Unless  the  input  signal  is  assumed  to  be 
white,  the  analysis  becomes  quite  complicated  due  to  the 
variation  of  the  variance  of  the  signal  and  therefore  of  the 
noise  sources  within  each  array.  Kancko  and  Liu  have  ob¬ 
tained  detailed  formulas  for  a  general  stochastic  model  of  the 
input  signal.  We  will  confine  attention  here  to  the  case  of 
white  input  signal,  where  the  signal  at  any  array  in  the  FFT 
is  also  white,  with  constant  variance  across  the  array. 

In  Fig.  12  a  typical  butterfly  computation  (only  top  half) 
is  indicated,  including  the  noise  sources  due  to  multiplication 
and  addition.  The  assumption  of  white  input  signal  implies 
that 

[Re  (A'*)]*  =  [Im  (,Ym)]-  =  ]  |A'mT  (68) 

and  application  of  our  floating-point  noise  model  as  in  Section 
III-D  yields  the  noise  source  variances 


<r„2  +  <t,52  =  <rfl2  +  <r,,:  =  <rfl2  =  a,,2  =  a,2|  -V„,|2  (69) 

O,2  =  a,, 2  =  a,2,.Ym[2.  (70) 
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The  variance  of  ihe  complex  noise  source  Um  =  is  then 


|trm|2  =  W|.Y,|*  (71) 

so  that  tlie  variance  of  the  noise  generated  in  computing  the 
(»i+l)th  array  is  4<r,2  times  the  variance  of  the  signal  in  the 
wth  array  If  the  input  (zeroth)  array  is  white  noise  with 
variance  a>2  then  the  noise  generated  in  the  (»/+l)th  array  is 
2m<rx2(4crt2) .  If  Onm*  is  the  output  noise  due  to  the  noise  gen¬ 
erated  in  the  (»i  +  l)tli  array,  then 
« 

<rom2  =  2-c»+'>2m(Tx2(4ff,2)  =  2Xa,W.  (72) 

Since  the  noise  ge  rated  in  each  array  is  assumed  to  he  inde¬ 
pendent,  the  total  output  noise  variance  or1  is 

or 2  =  2vA’o,!oI2.  (73) 

By  noting  that  the  output  signal  variance  is  related  to  the 
input  signal  variance  hy 


<r.v2  =  .\W  (74) 

the  result  follows: 

o 

air 

—  =  2°.2P.  (75) 

o\- 

A  further  result,  which  can  he  derived  from  our  model,  is 
an  expression  for  the  final  expected  output  noise-to-signal 
ratio  which  results  after  performing  an  FFT  and  an  inverse 
FFT  on  a  white  signal  ,v(«).  The  inverse  FFT  introduces  just 
as  much  roundolT  noise  as  the  FFT  itself,  and  thus  the  result¬ 
ing  output  noise-to-signal  ratio  is 

2 

—  ”  4ff, 2f  (76) 

<rx- 

or  just  double  the  result  in  (75). 

In  order  to  see  the  implications  of  (75)  or  (76)  in  terms  of 
register  length  requirements,  it  is  useful  to  express  these 
results  in  units  of  hits.  We  use 

(oir/ox^o,-)  bits  =  J  log2  ( 2v )  (77) 

to  represent  the  number  of  hits  hy  which  the  rms  noise-to- 
signal  ratio  increases  in  passing  through  a  floating  point  FFT. 
h'or  example,  for  i»  =  8  this  represents  2  hits  and  for  i>=  1 1  it 
represents  2.23  hits.  The  number  of  hits  of  rms  noise-to-signal 
ratio  increases  as  log>  (log*  N),  so  that  doubling  the  number 
of  points  in  the  FFT  produces  a  very  mild  increase  in  output 
noise,  significantly  less  than  the  half-bit-pcr-stage  increase 
for  fixed-point  computation.  In  fact,  to  obtain  a  half-hit  in¬ 
crease  in  the  result  above,  we  would  have  to  double  v,  or 
square  N. 

In  the  analysis  leading  to  (75),  we  have  not  considered  the 
fact  that  multiplications  hy  1  can  be  performed  noiselessly, 
h'or  a  specified  radix  2  algorithm,  such  as  the  decimation  in 
time  algorithm  shown  in  Fig.  10,  these  reduced  variances  for 
lT’  =  l  and  jean  he  included  in  the  model  to  obtain  a  slightly 
reduced  prediction  for  output  noise-to-signal  ratio.  However, 
for  reasonably  large  A-,  this  modified  noise  analysis  yields  only 
sliglifh  better  predictions  of  out  nit  noise  than  does  the 
simplified  analysis  above. 

A  consequence  of  our  analysis  leading  to  (75)  is  that  the 
output  noise  is  white  This  follows  from  the  fact  that  each 


Fig.  13.  (a)  Kxpcrimcma!  mid  theoretical  noise-to-signal  ratios  for  Hom¬ 

ing  |K)lnl  KKT.  (b)  ICxpcrimcntal  and  theoretical  noise-to-signal  ratios 
for  floating-point  KKT  and  inverse. 


array  of  noise  sources  is  white.  The  reduced  noise  source  vari¬ 
ance  for  lT  =  1  and  j  implies  that  for  some  arrays  there  will  he 
a  variation  of  noise  source  variance  over  the  array.  This  im¬ 
plies  a  slight  variation  of  output  noise  variance  over  the  out¬ 
put  array,  and  thus  our  modified  noise  analysis  will  only  pre¬ 
dict  an  average  noise  variance  over  the  output  array. 

The  results  discussed  above  have  been  verified  with  excel¬ 
lent  agreement  as  shown  in  Fig.  13(nj  and  (b).  To  obtain  this 
agreement,  however,  it  was  necessary  to  use  randomized 
rounding,  i.c.,  randomly  rounding  up  or  down  when  the  value 
of  mantissa  was  exactly  (1  2)2  *.  The  modified  theoretical 
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Fig.  t4.  F.X|>erimcnt:il  noi«i--l o-p»Kn:il  rallos  for  Homing  point  FKT  nnil 
I'i’T  -invL’Dv!  KKT;  iruncntion  used  insirad  of  rounding. 


curve  shown  was  obtained  by  taking  into  account  reduced 
noise  source  variances  for  If' «  1  and  If'  =  j.  Also  shown  arc 
experimental  results  for  nonrandomized  rounding.  These 
results  were  fitted  empirically  with  a  curve  of  the  form  «f*. 
but  this  quadratic  dependence  was  not  established  tlieo- 
r.  t’cally.  Noise-to-signal  ratios  were  also  measured  for  the 
ease  where  truncation  rather  than  rounding  was  used  in  the 
arithmetic;  the  results,  with  empirically  fitted  quadratic 
cun  es,  are  shown  in  Tig.  14. 

Our  analysis,  and  all  the  above  experiments,  applied  to  the 
case  of  white  signal.  Some  experimental  investigation  has 
been  carried  out  as  to  whether  the  predictions  are  valid  when 
the  signal  is  nonwhite.  Specifically,  the  noise  introduced  in 
computing  an  I’ FT  was  measured  for  sinusoidal  signals  of 
several  frequencies,  for  f-8.  9.  10,  and  11.  The  results,  aver¬ 
aged  over  the  input  frequencies  used,  were  within  15  percent 
of  those  predicted  by  (75).  In  these  ex |>cri incuts,  the  “ran¬ 
domized’’  rounding  procedure  was  used. 

E.  Effects  of  Coefficient  Quantization  in  the  EFT 

As  with  the  implementation  of  digital  filters,  the  imple¬ 
mentation  of  the  FFT  algorithm  requires  the  use  of  quantized 
coefficients.  While  a  completely  definitive  study  of  the  effects 
of  coefficient  quantization  in  the  FFT  remains  to  be  done,  two 
approaches  have  been  pursued  for  which  some  results  have 
been  obtained. 

Although  the  nature  of  coefficient  quantization  is  inher¬ 
ently  nonstatistical.  Weinstein  [24]  has  obtained  some  useful 
results  by  means  of  a  rough  statistical  analysis.  This  sta¬ 
tistical  analysis  corres|xmds  to  introducing  random  jitter  in 
the  coefficients  and  determining  tile  output  noise- to-signal 
ratio  due  to  this  noise.  While  the  detailed  effect  due  to  coeffi¬ 
cient  error  due  to  quantization  is  different  than  that  due  to  jit¬ 
ter,  it  is  reasonable  to  c.x|>ect  that  in  a  gross  sense  the  magni¬ 
tude  of  the  errors  is  comparable. 

To  develop  this  statistical  analysis,  we  let  F(k)  denote  the 
l)F1  of  a  sequence  f(n)  and  ffk)  the  result  of  transforming 
f(«)  with  a  radix  2  FFT  algorithm  with  jittered  coefficients. 


Then 


p(*)  =  I3/(»)IF"*  (78) 

n-0 

and 

Hh)  =  X)  /(»>)«.*.  (79) 

n-0 


Because  of  the  form  of  the  FFT  algorithm  each  element  11.* 
will  be  a  product  of  v-  log,  N  quantized  coefficients.  Thus 


¥ 

«-  -  II  (H  -  +  *.)  (80) 

»- 1 

where 

9 

II  IF-  -  IP*  (81) 

l-l 

with  b  bits  for  the  real  and  imaginary  parts  of  each  of  the 
coefficients,  excluding  sign,  |i(|  is  less  than  or  equal  to 
( v  2)2  If  we  assume  that  the  real  and  imaginary  parts  of 
the  jitter  in  the  coefficients  are  nncorrelated  and  uniformly 
distributed  lietween  plus  and  minus  ( 1/2)2 '*  then  <r*».  the 
variance  of  i,  is  a*’- 2~*/6.  The  error  in  the  computation  of 
the  DFT  can  be  expressed  as 

/•(A)  -  P(k)  -  E(k)  =  £ /(»)(».(  -  IF"*).  (82) 

Front  (80)  and  (81)  we  can  express  the  factor  (11,*-  IF"*)  as 

9  9 

(n.t  —  II  "*)  «=  23  8,  II  IF—  +  higher  order  terms.  (83) 

i-i  t-i 

If  we  neglect  higher  order  error  terms,  and  assume  that  5/ 
arc  mutually  uncorrelated  then  the  variance  of  (11.*-  IF"*)  is 
cqu.il  to  v(2~*/6).  Finally,  assuming  that  all  elements  Q.*  are 
uncorrelated  with  each  other  and  with  the  input  signal,  the 
output  error  variance  ok'  is 

**’-  „i_  £  | /(„)!».  (84) 

O  „«.o 

Since  from  ParsevaP*  relation 

El/Ml’-  4  El  '••(*)  I’ 

«-0  A  a-0 

«■  ntean-sqtiared  output  signal  (85) 

the  ratio  of  mean-squared  output  error  to  mean-squared  out¬ 
put  signal  is  thus 


ok 


Although  we  would  not  expect  (86)  to  predict  with  great 
accuracy  the  error  in  an  FFT  due  to  coefficient  quantization, 
it  is  helpful  as  a  rough  estimate  of  the  error  The  key  result  of 
(86),  which  we  woi  Id  like  to  test  ex|>erimcntally,  is  that  the 
err«r- to-signal  ratio  increases  very  mildly  with  .V.  living  pro- 
IHirnonal  to  K»logj  ,V,  so  that  doubling  <Y  produces  only  a 
slight  increase  in  the  error-to-signal  ratio. 


i  • 
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Fir.  IS.  Errors  due  to  coefficient  quantization  in  FFT. 

To  test  this  result,  experimental  measurements  on  errors 
due  to  coefficient  quantization  were  made.  In  each  run.  a 
sequence  /(»») — white  in  the  sense  that  all  2 N  real  numbers 
making  up  the  N- point  complex  sequence  were  mutually  un¬ 
correlated,  with  zero  means,  and  equal  variances — was  ob¬ 
tained  using  a  random  number  generator.  This  sequence  was 
transformed  twice,  once  using  a  36-bit  coefficient  table,  and 
once  using  a  coefficient  table  rounded  to  much  shorter  word 
length  (e.g.,  12  bits).  For  each  transform,  36-bit  accuracy  was 
used  in  the  arithmetic  to  make  the  effect  of  roundoff  error 
negligible.  The  results  were  subtracted,  squared,  averaged 
o\er  the  output  array,  and  divided  by  the  output  signal 
variance  ( A"  times  the  input  signal  variance)  to  obtain  an 
experimental  output  error-to-signal  ratio.  For  each  value  of 
A\  several  random  sequences  were  transformed  and  the  re¬ 
sults  were  averaged  to  obtain  statistically  convergent  esti¬ 
mates. 

These  results  are  displayed  in  Fig.  IS;  the  quantity  plotted 
is  ok1  l~*ort  where  oe1  is  the  mean-squared  output  signal  as 
defined  in  (85).  The  theoretical  curve  corresponding  to  (86)  is 
shown,  and  the  circles  represent  measured  output  error-to- 
signal  ratio  for  the  fixed-point  case.  We  note  that  the  "xperi- 
mental  results  generally  lie  below  the  theoretical  curve.  No 
experimental  result  differs  by  as  much  as  a  factor  of  two  from 
the  theoretical  result,  and  since  a  factor  of  two  in  oe'/of1  cor¬ 
responds  to  only  half-a-bit  difference  in  the  rnts  output  error, 
it  seemr  that  (86)  is  a  reasonably  accurate  estimate  of  the 
effect  of  coefficient  errors.  The  experimental  results  do  seem 
to  increase  essentially  linearly  with  v,  but  with  smaller  slope 
than  given  in  (86). 

In  the  above,  fixed-point  arithmetic  has  been  assumed. 
However,  since  a  block  floating-|mint  FFT  will  generally  use 
fixed-point  coefficients,  our  results  are  valid  for  the  block 
lloating-|Miint  case  also.  With  some  slight  modifications,  it  is 
possible  to  obtain  similar  results  for  the  floating-point  case. 
Except  for  a  constant  factor,  the  floating-  and  fixed-point 
results  are  the  same.  Experi mental  results  for  the  lloating- 
jsiint  case  are  represented  by  the  solid  dots  in  Fig.  15,  and 


are  observed  to  be  slightly  lower  than  the  results  foi  the  fixed- 
point  case. 

A  different  approach  to  the  characterization  of  FFT  coeffi¬ 
cient  quantization  has  been  taken  by  Tufts,  Hersey,  and 
Mosier  (29].  In  their  analysis  the  effect  of  coefficient  quantiza¬ 
tion  is  represented  in  terms  of  the  level  of  spurious  sidelobes 
introduced.  In  p.  'icular,  a  sequence  f(n)  =  «„(»  —  n„)  where 
«„(»)  denotes  a  unit  sample,  has  a  DFT  with  a  purely  sinus¬ 
oidal  real  and  imaginary  part,  i.e.. 

F(k)  =  IF"**  (87) 

and  the  inverse  DFT  of  F(k)  should,  of  course,  have  only  a 
single  nonzero  component.  Due  to  coefficient  quantization, 
however,  the  DFT  obtained  is 

F(k)  =  (88) 

and  since  the  real  and  imaginary  parts  are  not  exactly  sinus¬ 
oidal,  the  inverse  DFT,  with  exact  coefficients,  of  F(k)  will 
hnve  spurious  components.  For  each  of  the  set  of  N  sequences 

/*(»)  =  Uo(n  -  nk),  k  =  0,  1,  •  •  •  ,  X  —  1.  (89) 

Tufts  et  al.  compute  the  DFT  with  quantized  coefficients  fol¬ 
lowed  by  the  inverse  DFT  with  accurate  coefficients.  At  the 
output  of  the  inverse  DFT,  the  size  and  frequency  locations 
of  the  spurious  sidelohe  components  produced  due  to  thr 
quantized  coefficients  are  observed.  Since  any  function  f(n) 
can  lie  constructed  as  a  weighted  sum 

/(«)  =  Y  o*/*(")  (90) 

i-0 

the  spurious  sidclolies  produced  for  any  /(n)  can  in  principle 
be  determined  by  combining  the  responses  due  to  a  set  of  im¬ 
pulses.  But  currying  out  such  a  combination  is  not  practical 
for  arbitrary /(«).  Tufts  et  al.  have,  however,  tabulated  the 
worse  sidelohe  levels  encountered  for  any  }t(n)  as  a  function 
of  the  numlier  of  hits  retained  in  the  coefficients,  for  the  case  of 
a  64-point  FFT  and  sign-magnitude  representation  of  coeffi¬ 
cients. 

V.  Examples 

A.  Introduction 

In  the  preceding  sections  the  effects  of  arithmetic  roundoff 
have  been  analyzed  for  simple  (first-  and  second-order)  digital 
filters  and  the  FFT.  These  algorithms  arc  the  basic  building 
blocks  in  more  complicated  digital  processing  such  as  a  higher 
order  digital  filter  or  a  convolutional  filter  realized  via  the 
FFT.  Examples  will  lie  presented  in  this  section  to  indicate 
how  some  of  the  ideas  develop'd  above  can  lie  applied  to  ana¬ 
lyze  and  to  choose  the  most  advantageous  configuration  for 
such  systems.  The  first  two  examples  concern  the  realization 
of  higher  order  recursive  filters  and  Have  borrowed  from  the 
work  of  other  authors.  The  third  examp’  deals  with  an  FFT 
filter. 

h.  Fixed-Point  Digital  Filter  in  Cascade  and  Parallel  Form 

After  a  digital  filter  has  lieen  specified  in  terms  of  its  |k>Ics 
and  zeroes,  and  the  type  of  arithmetic  has  been  selected,  a 
choice  must  still  lie  made  among  the  various  possible  con¬ 
figurations  of  the  filter  which  will  differ  with  respect  to  the 
effects  of  roundoff  noise.  An  exhaustive  study  of  the  se  eetion 
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Fig.  16.  Second-order  section  with  [M)]cK  preceding  zeros.  Used  in  Jack- 
son’s  IP  parallel  form  with  fla"0.  Also  used  in  ID  direct  form. 


Fig.  I  7.  Second-order  section  with  zeros  preceding  {Males.  Used  in  2 P 
parallel  form  with  aj  — 0.  Also  used  in  lit  direct  iorm. 


of  filter  form  is  beyond  the  scope  of  this  paper,  blit  ait  excellent 
example  of  the  necessary  considerations  is  given  by  Jackson 
[30]  in  his  analysis  of  roundoff  noise  for  fixed-point  digital 
filters  realized  in  cascade  and  parallel  form 

Jackson  considers  two  parallel  form  realizations:  the  IP 
form  where  the  individual  second-order  sections  are  realized, 
as  shown  in  Fig.  16,  with  the  poles  preceding  the  zeros,  and 
the  2 P  form  where  zeros  precede  poles  in  the  individual  sec¬ 
tions.  as  shown  in  Fig.  17.  (Figs.  16  and  17  do  not  show  all  the 
scaling  coefficients  needed  to  prevent  overflow.)  His  analysis 
indicates  that  for  a  variety  of  scaling  criteria  (based  on  dif¬ 
ferent  Lp  norms*  of  the  input  signal)  and  for  various  measures 
of  the  output  noise  (such  as  its  total  power,  or  its  |x‘ak  s[x-c- 
tral  value),  the  output  signal-to-noise  ratios  of  the  two  forms 
are  very  close.  Generally,  a  very  slight  advantage  will  be 
gained  with  the  1 P  form. 

Comparison  of  the  two  parallel  forms  basically  reduces  to 
a  comparison  of  the  noise  properties  of  the  two  forms  of  sec¬ 
ond-order  seci  ons,  since  the  noises  from  the  second-order  sec¬ 
tions  are  simply  additive  in  the  output  of  the  parallel  form. 
Hence  the  discussion  above  applies  to  a  comparison  of  second- 
order  section  realizations.  Another  form  of  second-order  sec¬ 
tion  which  could  be  considered  is  a  coupled  form  as  shown  in 
Fig.  18.  For  the  case  a,  =  r  cos  8,  a-  =  r  sin  8,  this  filter  has 
poles  at  *  =  and  a  zero  at  i  =  r  cos  8.  The  coupled  form 
noise-to-signal  ratio  has  been  compared  (24)  to  ratios  for 
forms  essentially  the  same  as  those  in  Figs.  16  and  17  for  the 

•If 

*(•>  •  £  /(■>  eip  (  -jlw  — 

represents  the  Fourier  transform  of  a  signal  or  of  a  filler  impulse  resisinse 
/(a),  then  the  corresponding  l-p  norm  is 

where  u>,  denotes  sampling  frequency. 


®i 

y(n) 


•i 

Fig.  18.  Coupled  form  for  second-order  section. 

case  of  white-noise  input  and  an  absolute  overflow  constraint 
(through  the  type  of  analysis  given  in  Section  III-C).  The 
coupled  form  was  found  to  have  substantially  lower  noise-to- 
signal  ratio  for  filters  with  high  gain  and  low  resonant  fre¬ 
quencies.  For  l  =  1  —  r,  and  $^1,  the  results  vary  as 

<re2 

- ^  1/32  sin2  0  (forms  of  Figs.  16  and  17) 

<F,* 

<r,J 

- ~  1/4*  (coupleu  form).  (91) 

a,* 

The  implication  of  this  result,  together  with  the  somewhat 
reduced  coefficient  sensitivity  for  the  coupled  form,  is  that 
t‘  ' '»  form  may  lie  a  g<x*d  choice  in  some  situations,  despite  the 
fact  that  its  implementation  requires  four  multiplications 
instead  of  three  for  a  fx>le  pair  and  a  single  zero. 

As  state  a  above.  Jackson  found  not  much  difference  be¬ 
tween  the  noise  properties  of  the  1 P  and  2P  parallel  forms. 
However,  the  situation  more  interesting  in  the  case  of  the 
cascade  form.  Here  ik  finds  that  large  differences  are  possible 
between  the  roundoff  noise  outputs  of  the  1 D  (|x>les  before 
zeros  in  individual  sections)  and  ID  (zeros  liefore  jnilts  in 
each  section)  forms.  Also  the  ordering  of  the  sections  and  the 
pairing  of  |x>Ies  and  zeros  have  important  effect  on  the  output 
signa!-to-noise  ratio.  Jackson's  analyses  lead  to  several  rules 
of  thumb  for  selection  of  1 D  or  2D,  for  ordering  of  sections, 
and  for  pairing  of  |x>les  and  zeros. 

In  general,  the  choice  of  configuration  defends  on  which 
Lp  norm  of  the  scaled  transfer  functions  is  constrained  to  pi  *- 
vent  overflow  and  on  which  Lr  norm  of  the  output  noise  sjiec 
trum  is  used  as  a  measure  of  performance.  Two  Lp  constraints 
on  the  filter  are  of  particul  ir  intercs. :  the  />=  *  case,  where 
the  jieak  value  of  the  transfer  function  to  each  jxissible  over¬ 
flow  nod^  is  constrained;  and  the  p  —  2  case  where  the  rms 
transfer  function  to  each  ncxie  is  constrained.  The  choice 
*  is  just  slightly  less  stringent  than  the  absolute  overflow 
constraint  (7),  and  prevents  overflow  even  when  the  input  is 
a  narrow-band  signal  at  resonance  of  the  relevant  transfer 
function.  The  />  =  2  constraint  is  more  appropriate  for  prevent¬ 
ing  overflow  when  the  input  is  wide- band  in  nature.  Two  Lr 
norm^on  the  output  noise  spectrum  are  of  particular  interest: 
the  r=l  norm  which  measures  the  total  output  noise  power, 
and  the  norm  which  measures  the  |>eak  value  of  the 

spectrum  of  the  output  noise. 

With  regard  to  selection  of  \D  or  2D  forms,  Jackson's  rule 
of  thumb  says  to  select  1 D  when  2,  r*  «  and  2D  when 
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/>=oeIr=l;for^  =  2,  r  =  l  and  (or  p  =  »  ,  r  —  »  ,  either  form 
may  be  selected. 

The  choice  of  ordering  of  sections  also  de|H*nds  on  the 
norms  which  are  selected.  For  p  =  2,  r  *  «  ,  the  sections  should 
t>e  ordered  in  decreasing  |>eakedness,  where  peakedness  is 
defined  as  the  peak  gain  of  a  section  divided  by  its  rms  gain. 
For  p^  x  ,  r=  1,  the  sections  should  he  ordered  in  increasing 
peakedness.  For  p  =  2,  r=l  and  for  />=«,  r=<*.  the  choice 
of  ordering  depends  on  whether  form  2D  or  ID  is  chosen  for 
the  individual  sections.  Decreasing  peakedness  should  l>e 
chosen  with  form  ID,  and  increasing  peakedness  with  form 
ID. 

The  rule  for  pairing  of  poles  and  zeros  is  as  follows:  Let 
|  //.(«)  |  denote  the  magnitude  of  the  frequency  response  of 
the  »th  section,  and  M,  denote  the  maximum  over  u  of 
|  //»(&>)  | .  Then  the  pairing  should  be  chosen  such  that  the 
maximum  over  n  of  .1/.  is  minimized. 

The  above  rules  are  illustrated  by  Jackson  with  a  specific 
filter  example — a  sixth-order  Chebyshev  band  rejection  filter, 
and  the  results  are  in  accord  with  his  rules.  He  analyzes  the 
output  noise  of  this  filter  for  parallel  forms  f  P  and  IP  and  for 
all  orderings  of  cascade  forms  ID  and  2D  (with  proper  pole- 
zero  pairing).  Little  difference  is  seen  lietween  the  two  parallel 
forms.  For  p  =  2,  r  =  x  the  (teak  output  noise  spectrum  is  7-12 
dB  worse  for  the  ID  cascade  forms  than  for  the  ID  forms; 
while  for  p= x,  r=l,  the  output  noise  power  is  7-12  dB 
worse  for  ID  than  for  2D.  The  effects  of  |>< tit*  zero  pairing  and 
of  ordering  of  sections  also  follow  quite  well  the  rules  previ¬ 
ously  stated.  The  parallel  forms  turn  out  to  bo  slightly  supe¬ 
rior  to  the  liest  cascade  forms  with  resjiect  to  roundoff  noise. 

As  Jackson  indicates,  these  rules  of  thumb  have  certain 
qualifications  and  are  not  always  valid.  However,  they  have 
been  shown  to  be  helpful  in  a  variety  of  types  of  examples  [31  ]. 

C.  Choice  of  Form  for  Floating-Point  Digital  Filter 

By  means  of  an  example,  Liu  and  Kaneko  |fl]  have  com¬ 
pared  the  direct,  cascade,  and  parallel  form  realizations  of  a 
floating-point  digital  filter.  The  filter  selected  was  an  eighth- 
order  low-pass  elliptic  filter.  The  noise-to-signal  ratio  for  the 
parallel  form  was  about  20  dB  worse  than  for  the  direct  form, 
while  the  cascade  form  was  comparable  to  (alx>ut  1.5  dB 
worse  than)  the  parallel  form. 

Various  orderings  of  cascade  form  floating-point  filters 
have  not  l>een  studied  in  detail.  Probably  floating-point 
cascade  filters  are  not  too  sensitive  to  ordering  since  large 
variations  in  signal  level  from  stage  to  stage  can  lie  accom¬ 
modated  by  the  floating-point  exponent. 

A  comparison  of  the  noise-to-signal  ratio  properties  of 
floating-|N>int  second-order  sections  where  |siles  precede 
zeros  (Fig.  16)  and  where  zeros  |>rerede  |>oles  (Fig.  17)  indi¬ 
cates  that  at  least  for  white-noise  inputs  the  liehavior  of  the 
two  forms  is  essentially  identical.  For  a  high-gain  second  order 
section  of  low  resonant  frequency,  a  coupled  form  realization 
yields  some  noise-to-signal  ratio  advantage  over  both  of  these 
two  forms. 

D.  FFT  Filler 

The  results  of  our  roundoff  noise  analysis  for  lived  |M ti nt 
FFT  will  now  be  applied  to  obtain  an  expression  for  the  out¬ 
put  noise-to-signal  ratio  of  a  finite  impulse  response  digital 
filter,  implemented  by  means  of  the  FTT.  The  overflow  con¬ 
straints  of  this  ty|>e  of  filter  will  Ik-  accounted  for  in  the  analy¬ 


sis.  Attention  will  be  focussed  on  a  prototype  low-pass  filter 
with  256-point  impulse  response  and  a  cutoff  frequency  of  1/4 
the  half  sampling  frequency.  Rounded  arithmetic  will  be 
assumed. 

Let  us  examine  the  basic  steps  in  the  filtering  computation, 
tracing  the  buildup  of  noise  variance  as  we  proceed.  F'irst  the 
FFT  is  used  to  compute  the  DFT  of  a  section  of  input.  In  the 
implementation  of  a  filter  with  256-point  impulse  response,  it 
is  reasonable  to  compute  a  512-point  FFT,  where  the  input 
consists  of  256  data  samples  and  256  zeros.  Actually  512  real 
input  samples  would  be  treated  simultaneously,  by  placing 
sections  of  256  real  samples  in  both  the  real  and  imaginary 
parts  of  the  input  to  the  FFT.  To  guarantee  against  overflow, 
a  scaling  of  1/2  is  needed  at  each  stage  of  the  F'FT  yielding 
an  overail  attenuation  of  1/512.  The  samples  of  the  input  « 
sequence  must  be  less  than  unity  in  magnitude.  The  noise 
Ei(k)  at  the  output  of  this  first  DFT  has  variance 

®*i*  =  ]«.(*)  |*  =  |  2-**.  (92) 

This  noise  variance  is  small,  t>ecause  most  of  the  roundoff 
noise  has  been  shifted  off  by  the  attenuations.  However,  the 
scalings  have  also  caused  the  mean-squared  signal  to  decrease 
by  a  factor  of  1/512. 

Next  this  computed  transform  is  multiplied  by  a  sequence 
H(k)  representing  the  DFT  of  a  512-point  sequence  h(n)  con¬ 
sisting  of  the  filter  impulse  response  plus  256  zeros.  This  com¬ 
plex  multiplication  introduces  roundoff  noise  of  variance 
Tr*  5.  Assuming  that  we  have  chosen  '//(4)|  <1,  the  mean 
square  of  the  noise  E\  becomes  reduced  by 

I  //(*>  |*  =  B  (93) 

3 1 2  it  ,i 

a  ratio  of  the  filter  bandwidth  to  the  sampling  frequency. 

Thus  after  the  multiplication,  the  variance  of  ihe  total  noise 
Ei(k)  is 

2"26  5 

oK}  » - +  A*  -  2~ib.  (94) 

3  3 

This  noise  is  not  white,  but  has  a  component  whose  spectrum 
has  l**en  shaped  by  the  filter. 

For  the  example  under  consideration  1  4,  so 

3 

or.}  *»  -  2  **.  (95) 

4 

Note  that  or}  is  slightly  less  than  ok}  and  represents  or.ly  * 
about  a  bit  of  noise  However,  if  the  signal  spectrum  is  fiat,  * 
the  mean-squared  signal  will  also  Ik  reduced  somewhat  due 
to  the  multiplication  by  fl(k).  * 

Now  an  inverse  transform  is  computed  to  obtain  a  section 
of  output.  The  noise  variance  at  the  output  of  this  transform 
de|K»nds  on  how  many  scalings  are  necessary  in  the  inverse 
FFT.  In  order  to  detei  mine  how  many  scalings  are  necessary, 
a  lx>und  on  the  output  of  the  circular  convolution  is  required 
(32).  For  a  particular  filter,  such  a  U>und  can  Ik*  stated  as 

I  y(*) !  <  E  |  *0)  I  (96) 

M 

where  yin)  is  the  output  and  .1/  is  the  length  of  the  impulse 
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response.  The  prototype  filter  has  an  impulse  res|K>nse 

,  2  Tl  "  2ir kit 

h(n)  = - -  +  £(-!)*  cos - 

256 1_2  256 

2t(32)h  2  »(33)» 

+  0.7  cos -  0.225  cos - 


+  0.01995  cos  — — 


256 
2t(33)h- 
256  J 


256 


and 


(97) 


J 

A 


£  |  h(n)  |  =  3.12. 


(98) 


Hence,  only  two  scalings  (at  the  first  two  arrays)  are  neces¬ 
sary  in  the  inverse  transform.  Then,  in  propagating  through 
the  IFFT  (inverse  FFT),  the  variance  of  the  noise  £,(i) 
increases  by  a  factor  of  512/16.  (The  512  represents  the  gain 
of  the  inverse  DFT,  and  the  1/16  is  due  to  the  scalings.) 
The  variance  of  the  additional  output  noise  E,(k)  caused  by 
roundoff  in  the  IFFT  can  be  estimated  easily  via  the  method 
of  Section  1 V-C.  The  result  is 


on 


/512  51 2\  * 

M - + - )  +  o„'-  Y,  2*  =  (202)2-“. 

V  8  4  / 


(99) 


(100) 


The  total  mean-squared  output  noise  is 
512 

oh 2  =  —  oh.2  +  ok2  =  (226)2-“ 

16 

or  in  units  of  bits  of  .  ms  output  noise 

J  log,  (2 *W)  =  3.91  bits.  (101) 

The  mean-squared  output  signal  ran  be  estimated  if 
specific  statistics  are  assumed  for  the  input  signal  .v(»).  As  an 
example,  assume  that  *(n)  is  white  with  variance  <r,*  =  2  3. 
This  variance  goes  through  an  attenuation  of  1  512  in  the 
first  FFT,  an  attenuation  of  B  =  1  4  due  to  multiplication  by 
H(k),  and  a  gain  of  512  16  in  the  inverse  transform.  The 
mean-squared  output  signal  is  then 


and  the  output  noise- to-signal  ratio  is 

ok 2 

- =  (22  (MM))  2  “ 


(103) 


Assuming  an  input  noise-to-signal  ratio  (due*  to  A  I)  noise) 
of  (1  4)2  ^  the  noise-to-signal  ratio  has  worsened  by  a  factor 
of  alx>ut  55UO.  or 

\  log*  3500  =  6.1.5  bits  ( 1 04) 

ir  massing  through  the  FFT  filter. 
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